C
C *** LAST REVISED ON 11-JAN-1988 10:36:53.86
C *** SOURCE FILE: [DL.GRAPHICS.LONGLIB]MLIB2.FOR
C
C ************************************************************************
C
C MORE MASTER ROUTINES FOR THE LONGLIB GRAPHICS LIBRARY
C
C THESE ROUTINES ARE FORTRAN 77 COMPATIBLE WITH THE FOLLOWING EXCEPTIONS
C	1. LINES ARE INDENTED WITH TABS (^I).
C	2. TRAILING COMMENTS SEPARATE BY ! ARE USED.
C	3. A MACHINE DEPENDENT SCHEME USING %LOC() TO SIMULATE POINTERS
C	   IS USED IN VAX5D AND MVAX5D
C	4. CALLED ROUTINE NAMES AND COMMON BLOCKS MAY EXCEED 6 CHARACTERS.
C
C ************************************************************************
C
	SUBROUTINE LCNTR(V,NDX,NDY,NX,NY,NLEVELS,CVAL,NCH,
     1		NMIN,IFLAG,WORK,
     2		XAXL,YAXL,XTITL,NXT,YTITL,NYT,
     3		TITLE,NT,XMIN,XMAX,YMIN,YMAX,ICOL,ILINE)
C
C	CREATED BY D. LONG    MAY, 1986	AT JPL
C
C	MASTER ROUTINE TO PLOT A CONTOUR PLOT WITH LOG/LINEAR AXIS
C	THIS ROUTINE PLOTS CONTINUOUS CONTOURS AND CAN LABEL THEM
C	USING GCONTR
C
C	V	REAL DATA VALUE ARRAY DIMENSIONED V(NDX,NDY)
C	NDX,NDY	INT  DIMENSIONS OF V
C	NX	INT  NUMBER OF X VALUES
C	NY	INT  NUMBER OF Y VALUES
C	NLEVELS REAL NUMBER OF CONTOUR LEVELS
C			< 0 MIN(V) AND MAX(V) ARE USED WITH UNIFORM STEPS
C	CVAL	REAL ARRAY OF CONTOUR VALUES
C			IF NLEVELS < 0 THIS ARRAY IS FILLED WITH
C			THE VALUES USED FOR CONTOURS:
C			CVAL(J)=MIN(V)+(J-1)*(MAX(V)-MIN(V))/(J-1)
C	NCH	INT  LABELING OPTION FLAG
C			< 0 LABEL WITH CONTOUR VALUE (NCH IS NUMBER
C				OF DIGITS TO THE LEFT OF DECIMAL
C			= 0 NO LABELING OF CONTOURS
C			> 0 LABEL WITH LETTERS
C	NMIN	INT  MINIMUM NUMBER OF CELL FOR CONTOUR TO BE LABLED
C	IFLAG	INT
C		  (MAG) >10000 REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C			< 0  REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT)
C			= 0  CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C			> 0  SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C				NOTE: AXIS MAY BE LOG, BUT POINTS ARE PLOTTED 
C				      LINEARILY SPACED
C	  (ONE'S DIGIT) = 1 X AXIS LINEAR, Y AXIS LOGRITHMIC (BASE 10)
C	  (ONE'S DIGIT) = 2 X AXIS LOGRITHMIC,Y AXIS LINEAR
C	  (ONE'S DIGIT) = 3 X AXIS LOGRITHMIC,Y AXIS LOGRITHMIC
C	  (ONE'S DIGIT) = 4 X AXIS LINEAR, Y AXIS LINEAR
C	  (TEN'S DIGIT) = 0 NO GRID OR AXIS
C	  (TEN'S DIGIT) = 1 PLOT BOX WITH AXIS TICKS
C	  (TEN'S DIGIT) = 2 PLOT SOLID CARTESIAN GRID
C	  (TEN'S DIGIT) = 3 PLOT TICKED CARTESIAN GRID W/O BOX
C	  (TEN'S DIGIT) = 4 PLOT TICKED CARTESIAN GRID WITH BOX
C	  (TEN'S DIGIT) = 5 PLOT BOX WITH AXIS TICKS AND TICKED CARTESIAN GRID
C	  (TEN'S DIGIT) = 6 PLOT W/O BOX OR GRID BUT INCLUDE AXISES
C	  (TEN'S DIGIT) = 7 PLOT SOLID LOGARITHMIC GRID 
C	  (TEN'S DIGIT) = 8 PLOT DOTTED LOGARITHMIC GRID 
C	  (TEN'S DIGIT) = 9 PLOT TICKED LOGARITHMIC GRID 
C     (HUNDRED'S DIGIT) = 0 ASK WHICH SCREEN DEVICE TO USE
C		      <>0 SCREEN DEVICE CODE NUMBER
C	WORK	INT   WORKING ARRAY DIMENSIONED AT LEAST (2*NX*NY+1)/31
C	XAXL	REAL  X AXIS LENGTH IN INCHES
C			> 0  USE INPUT VALUE AXIS SCALING
C			< 0  USE SMOOTHED AXIS SCALING
C	YAXL	REAL  Y AXIS LENGTH IN INCHES
C			> 0  USE INPUT VALUE SCALING
C			< 0  USE SMOOTHED AXIS SCALING
C	XTITL	CHAR  X AXIS TITLE
C	NXT	INT   NUMBER OF CHARACTERS IN X AXIS TITLE STRING
C			< 0  TICKS ON INSIDE OF AXIS
C			= 0  NO AXIS
C			> 0  TICKS ON TITLE SIDE OF AXIS
C	YTITL	CHAR  Y AXIS TITLE
C	NYT	INT   NUMBER OF CHARACTERS IN Y AXIS TITLE STRING
C			< 0  TICKS ON INSIDE OF AXIS
C			= 0  NO AXIS
C			> 0  TICKS ON TITLE SIDE OF AXIS
C	TITLE	CHAR  PLOT TITLE
C	NT	INT   NUMBER OF CHARACTERS IN PLOT TITLE
C			= 0  NO TITLE
C			< 0  USE PEN COLOR ARRAY
C			IF ABS(NT)/100 > 0 THEN USE LINE TYPE ARRAY
C	XMIN	REAL  MINIMUM VALUE DISPLAYED ON X AXIS
C	XMAX	REAL  MAXIMUM VALUE DISPLAYED ON X AXIS
C	YMIN	REAL  MINIMUM VALUE DISPLAYED ON Y ARRAY
C	YMAX	REAL  MAXIMUM VALUE DISPLAYED ON Y ARRAY
C	ICOL	INT   ARRAY OF PEN COLORS DIMENSIONED ICOL(5+NLEVELS)
C			ICOL(1) = GRID COLOR
C			ICOL(2) = AXIS LINE COLOR
C			ICOL(3) = AXIS NUMBERS COLOR
C			ICOL(4) = AXIS TITLE COLOR
C			ICOL(5) = AXIS EXPONENT COLOR
C			ICOL(6) = TITLE COLOR (COLOR UPON RETURN)
C			ICOL(7) = CONTOUR LEVEL 1
C			ICOL(8) = CONTOUR LEVEL 2
C			   .        .      .   
C			ICOL(6+NLEVELS) = CONTOUR LEVEL NLEVELS (= AEND)
C	ILINE	INT  ARRAY OF PEN TYPES FOR CONTOURS
C
	INTEGER XTITL(1),YTITL(1),TITLE(1)
	REAL V(NDX,NDY),PY(2),CVAL(1)
	LOGICAL REPEAT
	INTEGER PL,ICOL(1),ILINE(1),IC(4)
	COMMON /CCNTRPLT/XM,DX,YM,DY
	DATA REPEAT/.FALSE./
	PL=3					!PRINTER DATA FILE= FOR003.DAT
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND !CLOSE PLOT PACKAGE
		REPEAT=.FALSE.
		RETURN
	ENDIF
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=-(JF/100)			!SCREEN DEVICE, CLEAR SCREEN
		CALL FRAME(PL,ILU,1.2,1.2,1.)	!INTIALIZE PLOT PACKAGE
	ELSE
		CALL CTERM(-1)			!RESTORE PLOT MODE
		JF=MOD(JF,100)
	ENDIF
	JF=MOD(JF,100)
	IF=JF
	IF (IF.GT.10) IF=MOD(IF,10)
	JF=JF-IF
	NGX=ANINT(ABS(XAXL))				!GRID SIZES
	NGY=ANINT(ABS(YAXL))
	PY(1)=XMIN
	PY(2)=XMAX
	CALL SCALE(PY,ABS(XAXL),2,1,1,XM2,DX2)
	IF (IF.EQ.2.OR.IF.EQ.3) THEN
		CALL SCALG(PY,ABS(XAXL),2,1,1,XM,DX)	!SMOOTH SCALE FACTORS
		IK=-1
		GDX=1./DX
		NGX=DX*ABS(XAXL)
		INX=-1
	ELSE
		CALL SCALE(PY,ABS(XAXL),2,1,1,XM,DX)	!SMOOTH SCALE FACTORS
		IK=0
		GDX=1.
		INX=1
	ENDIF
	IF (XAXL.GT.0.0) THEN				!INPUT SCALING
		XM2=XMIN
		DX2=(XMAX-XMIN)/ABS(XAXL)
		IF (IF.EQ.2.OR.IF.EQ.3) THEN
			XM=ALOG10(ABS(XMIN)+1.E-25)
			IF (XM.NE.AINT(XM)) XM=AINT(XM)-1.
			DX=AINT(ALOG10(ABS(XMAX)+1.E-25))+1.
			DX=(DX-XM)/ABS(XAXL)
			GDX=1./DX
			NGX=DX*ABS(XAXL)
		ELSE
			XM=XMIN
			DX=(XMAX-XMIN)/ABS(XAXL)
			GDX=1.
		ENDIF
	ENDIF
	PY(1)=YMIN
	PY(2)=YMAX
	CALL SCALE(PY,ABS(YAXL),2,1,1,YM2,DY2)
	IF (IF.EQ.1.OR.IF.EQ.3) THEN
		CALL SCALG(PY,ABS(YAXL),2,1,1,YM,DY)	!SMOOTH SCALE FACTORS
		IK=-1+IK
		GDY=1./DY
		NGY=DY*ABS(YAXL)
		INY=-1
	ELSE
		CALL SCALE(PY,ABS(YAXL),2,1,1,YM,DY)	!SMOOTH SCALE FACTORS
		IK=1+IK
		GDY=1.
		INY=1
	ENDIF
	IF (YAXL.GT.0.0) THEN				!INPUT SCALING
		YM2=YMIN
		DY2=(YMAX-YMIN)/ABS(YAXL)
		IF (IF.EQ.1.OR.IF.EQ.3) THEN
			YM=ALOG10(ABS(YMIN)+1.E-25)
			IF (YM.NE.AINT(YM)) YM=AINT(YM)-1.
			DY=AINT(ALOG10(ABS(YMAX)+1.E-25))+1.
			DY=(DY-YM)/ABS(YAXL)
			GDY=1./DY
			NGY=DY*ABS(YAXL)
		ELSE
			YM=YMIN
			DY=(YMAX-YMIN)/ABS(YAXL)
			GDY=1.
		ENDIF
	ENDIF
	IF (JF.EQ.0) GOTO 77		!NO AXIS
	IF (JF.EQ.60) GOTO 76		!NO GRID
	IF (JF.GT.60) THEN		!LOGRITHMIC GRID
		II=0
		IF (JF.EQ.70) II=1
		IF (JF.EQ.80) II=2
		CALL LGRID(0.,0.,GDX,GDY,INX*NGX,INY*NGY,II)
		GOTO 76
	ENDIF
	IF (JF.GE.20) THEN		!CARTESIAN GRID
		IF (JF.GE.30) NGX=-NGX
		IF (JF.GE.40) NGY=-NGY
		IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(1)),0.,0) !PEN COLOR
		CALL GRID(0.,0.,GDX,GDY,NGX,NGY)	!PRODUCE GRID
		IF (JF.EQ.50) JF=10
	ENDIF
76	NADD=0
	IF (NT.LT.0) THEN
		NADD=100000		!PEN COLOR
		IC(1)=ICOL(2)
		IC(2)=ICOL(3)
		IC(3)=ICOL(4)
		IC(4)=ICOL(5)
	ENDIF
	IX=IABS(NXT)
	XL=ABS(XAXL)*SIGN(1.,FLOAT(NXT))		!AXIS LENGTH WITH TICK INFOR
	IF (NXT.NE.0) THEN
		IF (IF.EQ.2.OR.IF.EQ.3) THEN
		CALL LGAXS(0.,0.,XTITL,-IX-NADD,XL,0.,XM,DX,IC) !X-AXIS
		IF (JF.EQ.10) CALL LGAXS(0.,ABS(YAXL),XTITL,IX+100+NADD,
     $			XL,0.,XM,DX,IC) !X-AXIS
		ELSE
		CALL AXIS(0.,0.,XTITL,-IX-NADD,XL,0.,XM,DX,N1,N2,IC)  !X-AXIS
		IF (JF.EQ.10) CALL AXIS(0.,ABS(YAXL),XTITL,IX+100+NADD,
     $			XL,0.,XM,DX,N1,N2,IC)  !X-AXIS
		ENDIF
	ENDIF
	IY=IABS(NYT)+1000
	YL=ABS(YAXL)*SIGN(1.,FLOAT(NYT))
	IF (NYT.NE.0) THEN
		IF (IF.EQ.1.OR.IF.EQ.3) THEN
		CALL LGAXS(0.,0.,YTITL,IY+NADD,YL,90.,YM,DY,IC) !Y-AXIS
		IF (JF.EQ.10) CALL LGAXS(ABS(XAXL),0.,YTITL,-IY-100-NADD,
     $			YL,90.,YM,DY,IC) !Y-AXIS
		ELSE
		CALL AXIS(0.,0.,YTITL,IY+NADD,YL,90.,YM,DY,N1,N2,IC)  !Y-AXIS
		IF (JF.EQ.10) CALL AXIS(ABS(XAXL),0.,YTITL,-IY-100-NADD,
     $			YL,90.,YM,DY,N1,N2,IC)	!Y-AXIS
		ENDIF
	ENDIF
C
77	CONTINUE
	NL=IABS(NLEVELS)
	IF (NL.EQ.0) NL=2
	AE=1.E20
	IF (NLEVELS.LE.0) THEN
		AE=V(1,1)
		AS=V(1,1)
		DO 65 IX=1,IABS(NX)
			DO 65 IY=1,NY
				AE=MAX(AE,V(IX,IY))
				AS=MIN(AS,V(IX,IY))
65		CONTINUE
		DO 66 IX=1,NL
			CVAL(IX)=AS+(IX-1)*(AE-AS)/FLOAT(NL-1)
66		CONTINUE
	ENDIF
	DXX=XAXL/FLOAT(IABS(NX-1))
	DYY=YAXL/FLOAT(IABS(NY-1))
C
	ICL=0
	IF (IABS(NT)/100.GT.0) ICL=2
C
C	CALL CONTOURING ROUTINE ONCE FOR EACH CONTOUR LEVEL TO MINIMIZE
C	STORAGE SPACE AND PEN CHANGES
C
	DO 100 IX=1,NL
		CV=CVAL(IX)
		IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(IX+6)),0.,0)	! PEN COLOR
		IF (ICL.NE.0) IC(1)=ILINE(IX)			! LINE TYPE
		CALL GCONTR(V,NDX,NDY,NX,NY,DXX,DYY,CV,-IX,AE,WORK,
     $			NCH,0.15,NMIN,ICL,IC,IC(1))
100	CONTINUE
C
C	PLOT TITLE AND FINISH UP ROUTINE
C
	IF (IABS(NT)/100.GT.0) CALL NEWPEN(0)		!RESET LINE TYPE
	CS=.21						!TITLE CHARACTER SIZE
	IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(6)),0.,0)	!PEN COLOR
	IF (NT.NE.0.AND.JF.NE.0) THEN			!TITLE
		XP=0.0
		YP=0.0
		CALL SYMBOL(XP,YP,CS,TITLE,
     $			0.,MOD(IABS(NT),100),-3)	!TITLE LENGTH
		XP=ABS(XAXL/2.)-XP/2.
		IF (XP.LT.0.) XP=0.
		CALL SYMBOL(XP,ABS(YAXL)+.15,CS,TITLE,
     $			0.,MOD(IABS(NT),100),-1)	!PLOT TITLE
	ENDIF
	CALL PLOT(0.,0.,3)
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)				!ASK IF CLEAR SCREN
		CALL PLOTND				!CLOSE PLOT PACKAGE
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)				!RETURN TO TEXT MODE
	ENDIF
	RETURN
	END
C
C *******************************************************************
C
	SUBROUTINE CNTLN(X,Y,Z,N,XH,YH,IFLAG,NCNTRS,CLIST,IAXIS,
     $		XT,NXT,TX,TSX,FSX,YT,NYT,TY,TSY,FSY,
     $		T,NT,XAST,XAEN,YAST,YAEN,ICOL,ILINE)
C
C	CREATED BY D. LONG     FEB, 1986	AT JPL
C
C	ROUTINE TO PLOT DATA CONTOURS SPECIFIED AS (X,Y,Z) TRIPLES 
C
C	X,Y,Z	(R) ARRAYS CONTAINING DATA POINTS
C	N	(I) NUMBER OF DATA POINTS
C	XH,YH   (R) LENGTH OF EACH AXIS (INCHES)
C	IFLAG	(I) PLOT FLAG
C		>10000	REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0	REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT)
C		= 0	CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0	SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C 	(ONE'S DIGIT) = 1 DO NOT USE COLOR ARRAY OR LINE TYPE ARRAY
C		      = 2 USE COLOR ARRAY BUT NOT LINE TYPE ARRAY
C		      = 3 DO NOT USE COLOR ARRAY BUT USE LINE TYPE ARRAY
C		      = 4 USE COLOR ARRAY AND LINE TYPE ARRAY
C 	(TEN'S DIGIT) = 0 JUST X,Y AXES
C		      = 1 X AND Y AXES IN A BOX
C 	(HUNDRED'S)   = 0 ASK WHICH SCREEN DEVICE TO USE
C		      <>0 SCREEN DEVICE CODE NUMBER
C	NCNTRS	(I) NUMBER OF CONTOUR LINES TO BE DRAWN
C		    (SELF COMPUTING IF NCNTRS.LE.0)
C	CLIST   (R) ARRAY OF CONSTANT CONTOUR VALUES IF NCNTRS > 0
C		NOTE: IF NCNTRS <= 0, THEN CLIST(1) = BASE VALUE,
C		AND CLIST(2) = INCREMENT VALUE (DELTA)
C 	IAXIS	(I) AXIS OPTION FLAG
C			< 0 DO NOT PLOT AXES
C 	(ONE'S DIGIT)	= 1 PLOT Y AXIS USING COMPUTED SCALE
C			= 2 PLOT Y AXIS USING COMPUTED, SMOOTHED SCALE
C			= 3 PLOT Y AXIS USING INPUT SCALE--VARIABLES REQUIRED
C			= 4 PLOT Y AXIS USING INPUT, SMOOTHED SCALE--VARIABLES REQUIRED
C 	(TEN'S DIGIT)	= 1 PLOT X AXIS USING COMPUTED SCALE
C			= 2 PLOT X AXIS USING COMPUTED, SMOOTHED SCALE
C			= 3 PLOT X AXIS USING INPUT SCALE--VARIABLES REQUIRED
C			= 4 PLOT X AXIS USING INPUT, SMOOTHED SCALE--VARIABLES REQUIRED
C	(HUNDRED'S)	= 0 NORMAL
C			= 1 TRIANGULATION ONLY PLOTTED
C 	XT,YT	(B) STRINGS FOR AXIS TITLES
C 	NXT,NYT	(I) LENGTH OF AXIS TITLES
C		     IF ZERO THEN AXIS NOT PLOTTED
C		     IF LESS THAN ZERO THEN AXIS TICKS ARE OPPOSITE AXIS TITLE
C FOLLOWING ONLY REQUIRED IF IAXIS > 0
C	TX,TY	(R) NUMBER OF TICK MARKS (SEE AXIS3)
C	TSX,TSY	(R) SIZE OF TITLE AND NUMBERS OF AXIS
C			  IF LESS THAN ZERO DO NOT AUTO-SCALE BY (x10^POWER)
C	FSX,FSY	(R) NUMBER FORMAT OF AXISES (SEE AXIS3)
C	T	(B) PLOT TITLE
C	NT	(I) NUMBER OF CHARACTERS IN TITLE
C			IF NT=0 NO TITLE PLOTTED
C 	XAST,YAST(R) AXIS START VALUES
C 	XAEN,YAEN(R) AXIS END VALUES
C	ICOL	(I) COLOR INDEX TABLE (REQUIRED IF 1'S DIGIT OF IFLAG=2
C			ICOL(1) AXIS LINE COLOR
C			ICOL(2) AXIS NUMERIC LABEL COLOR
C			ICOL(3) AXIS LABEL COLOR
C			ICOL(4) AXIS EXPONENT COLOR
C			ICOL(5) TITLE COLOR
C			ICOL(6) SURFACE COLOR
C	ILINE	(I) LINE TYPE ARRAY
C
C	EXTERNAL CALLS:
C
C		INTERPC
C		TRIANGC
C		  MIDC (FUNCTION)
C		CNTOUR
C		  CNTCRV
C ---------------------------------------------------------------------
C
	DIMENSION X(1),Y(1),Z(1)
	DIMENSION CLIST(2),ICOL(1),ILINE(1),IC(4)
	INTEGER XT(1),YT(1)
	LOGICAL REPEAT/.FALSE./
C
	DIMENSION PAS(2)
	DIMENSION IT(1000,3),IP(500),IB(500)
	DATA NZZ/1000/,NW/500/,NZ/1494/
	DIMENSION  IE(1494,2),ITE(1494,4),XI(1494),ETA(1494),
     *		LAMBDA(1494),IBE(1494)
C
	COMMON /CCNTRPLT/XMIN,DXR,YMIN,DYR
C
	IF (N.GT.NW) THEN
		CALL CTERM(1)
		WRITE(*,3010)
3010		FORMAT(' *** WARNING: INSUFFICIENT WORKING SPACE')
		CALL CTERM(-1)
	ENDIF
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	XLEN=ABS(XH)
	YLEN=ABS(YH)
	IYAXIS=MOD(IABS(IAXIS),10)
	IXAXIS=MOD(IABS(IAXIS),100)/10
	ITRIANG=MOD(IABS(IAXIS),1000)/100
	IF (IXAXIS.GT.2) THEN
		XMAX=XAEN
		XMIN=XAST
	ELSE
		XMAX=X(1)
		XMIN=X(1)
		DO 16 I=1,N
			XMAX=AMAX1(XMAX,X(I))
			XMIN=AMIN1(XMIN,X(I))
16		CONTINUE
	ENDIF
	IF (IXAXIS.EQ.2.OR.IXAXIS.EQ.4) THEN
		PAS(1)=XMAX
		PAS(2)=XMIN
		CALL SCALE(PAS,XLEN,2,1,1,XMIN,DAA)
		XMAX=YLEN*DAA+XMIN
	ENDIF
C
	IF (IYAXIS.GT.2) THEN
		YMAX=YAEN
		YMIN=YAST
	ELSE
		YMAX=Y(1)
		YMIN=Y(1)
		DO 17 I=1,N
			YMAX=AMAX1(YMAX,Y(I))
			YMIN=AMIN1(YMIN,Y(I))
17		CONTINUE
	ENDIF
	IF (IYAXIS.EQ.2.OR.IYAXIS.EQ.4) THEN
		PAS(1)=YMAX
		PAS(2)=YMIN
		CALL SCALE(PAS,YLEN,2,1,1,YMIN,DAA)
		YMAX=YLEN*DAA+YMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=MOD(JF,10000)/100
		CALL FRAME(3,ILU,1.5,0.65,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL
 	ENDIF
	JFF=MOD(JF,100)/10
	JF=MOD(JF,10)
C
	IF (IAXIS.GT.0) THEN				!PLOT AXIS LABELS
		NADD=0
		IF (JF.EQ.2.OR.JF.EQ.4) THEN
			NADD=100000
			IC(1)=ICOL(1)
			IC(2)=ICOL(2)
			IC(3)=ICOL(3)
			IC(4)=ICOL(4)
		ENDIF
		IF (NXT.NE.0) THEN 	!PLOT X AXIS
			CALL AXIS3(0.,0.,XT,-NADD-IABS(NXT),
     $				XLEN*ISIGN(1,NXT),
     $				0.,XMIN,XMAX,TX,TSX,FSX,IC)
			IF (JFF.EQ.1) CALL AXIS3(XLEN,0.,XT,
     $				NADD+101,XLEN*ISIGN(1,NXT),0.,
     $				XMIN,XMAX,TX,TSX,FSX,IC)
		ENDIF
		IF (NYT.NE.0) THEN 	!PLOT Y AXIS
			CALL AXIS3(0.,0.,YT,IABS(NYT)+1000+NADD,
     $				YLEN*ISIGN(1,NYT),
     $				90.,YMIN,YMAX,TY,TSY,FSY,IC)
			IF (JFF.EQ.1) CALL AXIS3(0.,0.,YT,101+NADD,
     $				YLEN*ISIGN(1,NYT),
     $				90.,YMIN,YMAX,TY,TSY,FSY,IC)
		ENDIF
	ENDIF
C
C	ADD TITLE
C
	IF (JF.EQ.2.OR.JF.EQ.4) CALL PLOT(FLOAT(ICOL(5)),0.,0) ! COLOR
	IF (NT.GT.0) THEN
		X1=0.0
		Y1=0.0
		CALL SYMBOL(X1,Y1,0.25,T,0.,NT,-3)	! TITLE LENGTH
		X1=(XLEN-X1)/2.
		CALL SYMBOL(X1,YLEN+.1,0.25,T,0.,NT,-1)	! PLOT TITLE
	ENDIF
C
C	BEGIN SURFACE PLOTTING
C
	DX=1.0
	DY=1.0
	IF (XMAX-XMIN.NE.0.0) DX=XLEN/(XMAX-XMIN)
	IF (YMAX-YMIN.NE.0.0) DY=YLEN/(YMAX-YMIN)
C
C CALL ROUTINE TRAING TO TRIANGULATE X-Y DATA POINTS
C
	CALL TRIANGC(X,Y,N,IT,NZZ,M,IP,IB,NW,LEDGES,NZ,IE,IBE,ITE)
C
C	PLOT TRIANGULATED SURFACE IF DESIRED
C
	DXR=1./DX
	DYR=1./DY
	IF (ITRIANG.EQ.1) THEN
		DO 35 I=1,M
			IPP=3
			DO 33 J=1,3
				L=IT(I,J)
				X1=DX*(X(L)-XMIN)
				Y1=DY*(Y(L)-YMIN)
				CALL PLOT(X1,Y1,IPP)
				IPP=2
33			CONTINUE
			L=IT(I,1)
			X1=DX*(X(L)-XMIN)
			Y1=DY*(Y(L)-YMIN)
			CALL PLOT(X1,Y1,2)
35		CONTINUE
		GOTO 300
	ENDIF
C
C DETERMINE THE RANGE OF THE Z DATA UNDER CONSIDERATION
C
  120	ZMIN = Z(1)
	ZMAX = ZMIN
	DO 50 K=2,N
		ZMIN = AMIN1(ZMIN,Z(K))
		ZMAX = AMAX1(ZMAX,Z(K))
   50	CONTINUE
	K=0
C
C HAS A CONTOUR LIST BEEN GIVEN? . .
C
	K=0
	FK = -1.0
  200	IF (NCNTRS.GT.0) GOTO 180
C
C DETERMINE (NEXT) CONTOUR CONSTANT VALUE
C
  210	FK = FK+1.0
	ZCON = FK*CLIST(2) + CLIST(1)
	IF (ZCON.GT.ZMIN.AND.ZCON.LT.ZMAX) GOTO 150
	IF (ZCON.GT.ZMAX.AND.CLIST(2).GE.0.) GOTO 300
	IF (ZCON.LT.ZMIN.AND.CLIST(2).LE.0.) GOTO 300
	GOTO 210
C
  180	K = K+1
	IF (K.GT.NCNTRS) GOTO 300
	ZCON = CLIST(K)
	IF (ZCON.LT.ZMIN.OR.ZCON.GT.ZMAX) GOTO 200
C
C CALL ROUTINE INTERPC TO
C INTERPOLATE FOR CONTOUR LINE DATA POINTS
C
  150	CALL INTERPC(X,Y,Z,N,ZCON,LEDGES,IE,LAMBDA,XI,ETA,J,NZ)
C
C ANY DATA POINTS FOUND? . .
C CALL ROUTINE CNTOUR TO SORT THE INTERPOLATED POINTS
C ON THE CONTOUR LINE AND DRAW IT
C
	IF (JF.EQ.2.OR.JF.EQ.4) CALL PLOT(FLOAT(ICOL(K+5)),0.,0)! COLOR
	IF (JF.EQ.3.OR.JF.EQ.4) CALL NEWPEN(ILINE(K))		! LINE TYPE
	IF (J.NE.0) CALL CNTOUR(ZCON,XI,ETA,LAMBDA,J,IBE,ITE,NZ,
     $		DX,XMIN,DY,YMIN)
	GOTO 200
C
C	FINISH ROUTINE
C
  300	CALL PLOT(0.,0.,3)			!PEN UP
	IF (JF.EQ.3.OR.JF.EQ.4) CALL NEWPEN(0)	!RESET LINE TYPE
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL TEXT MODE
	ENDIF
	RETURN
	END
C
C
	SUBROUTINE INTERPC (X,Y,U,N,ZCON,LEDGES,IE,LAMBDA,XI,ETA,J,NZ)
C
C	COSMIC CONTOUR ROUTINES ARC-11441
C	MODIFIED  D.G. LONG, JPL  FEB 1986
C ------------------------------------------------------------------
C
C ROUTINE INTERPC IS GIVEN A CONSTANT U VALUE (BIGU) FOR WHICH
C THE CONTOUR LINE IS TO BE DRAWN.  CHECK ALL GIVEN TRIANGLE EDGES,
C (ARRAY IE) AND CHECK THE VALUES OF U AT THE ENDPOINTS.
C INTERPOLATE FOR ALL POSSIBLE VALUES ON THE TRIANGLE EDGES.
C USE A LINEAR INTERPOLATION SCHEME
C
C X,Y,U     = DEPENDENT AND INDEPENDENT VALUES FOR
C             THE RELATION U=F(X,Y)  (INPUT)
C ZCON      = CONSTANT VALUE OF Z FOR WHICH INTERPOLATION
C             IS REQUIRED  (INPUT)
C LEDGES    = NO. OF EDGES IN THE TRIANGULATION  (INPUT)
C IE        = EDGE ENDPOINT INDICES FROM TRIANGULATION (INPUT)
C LAMBDA    = INDEX OF EDGES FOR INTERPOLATED POINTS  (RETURN)
C XI        = LIST OF X-COORDINATES OF INTERPOLATED POINTS
C ETA       = LIST OF Y-COORDINATES OF INTERPOLATED POINTS
C J         = NUMBER OF VALUES IN XI, ETA LISTS
C             (XI, ETA AND J ARE RETURNED)
C NZ	    = DIMENSION OF IE,XI,ETA,LAMBDA
C     ------------------------------------------------------------------
C
	DIMENSION  X(N),Y(N),U(N)
	DIMENSION  IE(NZ,2),XI(NZ),ETA(NZ),LAMBDA(NZ)
C
	J = 0
C
	DO 1 LCNT=1,LEDGES
C
C (A)
C DETERMINE X,Y,Z FOR THE ENDPOINTS OF THE NEXTEDGE - ORDER THEM
C
	I1 = IE(LCNT,1)
	I2 = IE(LCNT,2)
	X1 = X(I1)
	X2 = X(I2)
	Y1 = Y(I1)
	Y2 = Y(I2)
	U1 = U(I1)
	U2 = U(I2)
C
C (B)
C FUNCTION VALUES EQUAL AT ENDPOINTS OR
C CONSTANT ZC NOT BETWEEN THEM? . .
C
	IF (U1.EQ.U2) GOTO 1
	IF (U1.LT.U2) GOTO 100
	TEMP = U2
	U2 = U1
	U1 = TEMP
	TEMP = X2
	X2 = X1
	X1 = TEMP
	TEMP = Y2
	Y2 = Y1
	Y1 = TEMP
  100	IF (ZCON.LT.U1.OR.U2.LT.ZCON) GOTO 1
	IF (U2.EQ.ZCON) U2 = 1.000001 * ZCON
	J = J+1
C
C (E,F)
C LINEAR INTERPOLATION IS REQUIRED
C FOR THIS EDGE OVER THE Z-SURFACE
C
  101	T1 = (U2-ZCON)/(U2-U1)
	T2 = (ZCON-U1)/(U2-U1)
	XI(J) = T1*X1+T2*X2
	ETA(J)= T1*Y1+T2*Y2
  200	LAMBDA(J) = LCNT
    1	CONTINUE
	RETURN
	END
C
	SUBROUTINE CNTOUR (ZCON,XI,ETA,LAMBDA,J,IBE,ITE,NZ,DX,XM,DY,YM)
C
C	COSMIC CONTOUR ROUTINES ARC-11441
C	MODIFIED  D.G. LONG, JPL  FEB 1986
C ------------------------------------------------------------------
C
C A SET OF J INTERPOLATED POINTS FOR Z=ZCON  (XI(I),ETA(I) ON EDGE
C LAMBDA(I) FOR I=1,J), THE CONTOUR LINES MUST NOW BE DRAWN.  THERE
C MAY BE SEVERAL LINES, EITHER OPEN OR CLOSED CONTOURS.  THIS
C ALGORITHM WILL USE THE TRIANGULATION RELATIONSHIPS TO SORT OUT
C EACH LINE IN ORDER.  AS EACH CONTOUR LINE IS ESTABLISHED, USER
C SUPPLIED PROGRAM CNTCRV IS CALLED TO OUTPUT IT TO THE GRAPHICS
C DEVICE BEING USED.
C
C ARGUMENTS (ALL ARE INPUTS) -
C  ZCON    = CONSTANT VALUE OF Z UNDER CONSIDERATION
C  XI(J)   = ARRAY OF X COORDIANTES OF INTERPOLATED POINTS
C  ETA(J)  = ARRAY OF Y COORDIANTES OF INTERPOLATED POINTS
C  LAMBDA(J) = ARRAY OF EDGE NUMBERS FOR J-TH INTERPOLATED POINT
C  J       = NUMBER OF POINTS IN THE LIST OF INTERPOLATED POINTS
C  IBE     = THE LIST OF BOUNDARY EDGES TAKEN FROM THE TRIANGULATI
C  ITE     = LINKED LIST OF INDICES OF ADJACENT EDGES PROVIDED
C  		BY THE TRIANGULATION PROCEDURE.
C  NZ	   = DIMENSION OF XI,ETA,LAMBDA,IBE,ITE
C  DX,XM   = PLOTTING SCALE FACTOR AND MINIMUM VALUE FOR X DIRECTION
C  DY,YM   = PLOTTING SCALE FACTOR AND MINIMUM VALUE FOR X DIRECTION
C
C  EXTERNAL CALLS:
C   (CNTCRV)/(PLOT)
C
C ------------------------------------------------------------------
C
	DIMENSION  XI(NZ),ETA(NZ),LAMBDA(NZ),IBE(NZ),ITE(NZ,4)
CCC	DIMENSION XX(1495),YY(1495)
C
C (A)
C INITIALIZE LOCAL VARIABLES
C
	IF (J.EQ.0) RETURN
   10	J1 = 0
C
C (B,C)
C SEARCH THE LIST OF EDGES FOR A BOUNDARY EDGE (BE(I)=1)
C
    1	J1 = J1+1
	L1 = LAMBDA(J1)
	IF (IBE(L1).EQ.1) GOTO 2
	IF (J1.LT.J) GOTO 1
	GOTO 11
C
C SEARCH FOR A BOUNDARY EDGE AND PUT IT AT THE TOP OF THE LIST.
C
C (D)
C PUT THIS INTERPOLATED POINT AT THE TOP OF THE
C LIST FOR THIS CONTOUR,  SET J1
C
    2	IF (J1.EQ.J) GOTO 3
	XI(J+1) = XI(J1)
	ETA(J+1) = ETA(J1)
	LAMBDA(J+1) = LAMBDA(J1)
	DO 101 JCNT = J1,J
		XI(JCNT) = XI(JCNT+1)
		ETA(JCNT) = ETA(JCNT+1)
  101		LAMBDA(JCNT) = LAMBDA(JCNT+1)
C
C (E)
C SEARCH THE REMAINING POINTS FOR AN ADJACENT (COMMON) EDGE
C
    3	J1BIG = J
	LCNT = L1
    6	J1BIG = J1BIG-1
	J1 = 0
    5	J1 = J1+1
	L1 = LAMBDA(J1)
	DO 102 I=1,4
	IF (L1.EQ.ITE(LCNT,I)) GOTO 4
  102	CONTINUE
C (F)
C ERROR - THERE IS NO NEXT POINT.
	IF (J1.LT.J1BIG) GOTO 5
	GOTO 800
C
C (G)
C PUT THIS POINT AT THE TOP OF THE LIST.
C CONTINUE IF ITS NOT A BOUNDARY EDGE.
C
    4	XI(J+1) = XI(J1)
	ETA(J+1) = ETA(J1)
	LAMBDA(J+1) = LAMBDA(J1)
	DO 103 JCNT = J1,J
		XI(JCNT) = XI(JCNT+1)
		ETA(JCNT) = ETA(JCNT+1)
  103		LAMBDA(JCNT) = LAMBDA(JCNT+1)
	LCNT = L1
	IF (IBE(L1).NE.1) GOTO 6
C
C (H)
C DRAW THE OPEN CONTOUR LINE THROUGH THE POINTS
C XI(J1),ETA(J1) ...... XI(J1+1),ETA(J1+1) ...... XI(J),ETA(J)
C THEN RESET J AND CONTINUE
C ----------------------------------------------------------------
C
	NPOINT = J-J1BIG+1
	IF (NPOINT.LE.1) GOTO 300
	IP = 3
	DO 290 I=1,NPOINT
		X1=DX*(XI(I+J1BIG-1)-XM)
		Y1=DY*(ETA(I+J1BIG-1)-YM)
		CALL PLOT(X1,Y1,IP)
		IP = 2
290	CONTINUE
CC	CALL CNTCRV (XI(J1BIG),ETA(J1BIG),NPOINT,ZCON)
C
C ----------------------------------------------------------------
C
  300	J=J1BIG - 1
C
C (I)
C ARE THERE ANY MORE POINTS LEFT? . .
C
	IF (J) 800,800,10
C
C (J)
C NOW DRAW INTERNAL LINES (CLOSED CONTOURS THAT DO NOT START
C OR STOP AT BOUNDARY EDGES).  THE POINT AT J1BIG=J IN
C THE LIST IS CHOSEN TO START THE CONTOUR.
C
   11	J1BIG = J+1
	LCNT = LAMBDA(J)
C
C (K,M,P)
C FIND THE NEXT POINT FOR THIS CONTOUR (ON AN EDGE WITH A COMMON
C END POINT).  PUT IT AT THE TOP OF THE LIST, AND REPEAT UNTIL
C NO MORE COMMON EDGES REMAIN FOR THIS LINE.
C
   16	J1BIG = J1BIG-1
	J1 = 0
	IF (J1BIG.GT.J) J1=1
   15	J1 = J1+1
	L1 = LAMBDA(J1)
	DO 104 I=1,4
		IF (L1.EQ.ITE(LCNT,I)) GOTO 14
  104	CONTINUE
	IF (J1.LT.J1BIG) GOTO 15
C (L)
C OTHERWISE, NO ADJACENT EDGE WAS FOUND.
C THIS CONTOUR LINE IS COMPLETE, GO DRAW IT.
	GOTO 17
C
   14	XI(J+1) = XI(J1)
	ETA(J+1) = ETA(J1)
	LAMBDA(J+1) = LAMBDA(J1)
	DO 105 JCNT = J1,J
		XI(JCNT) = XI(JCNT+1)
		ETA(JCNT) = ETA(JCNT+1)
  105		LAMBDA(JCNT) = LAMBDA(JCNT+1)
	LCNT = L1
	IF (J1BIG.NE.1) GOTO 16
C
C (O)
C DRAW THE CLOSED CONTOUR LINE, THE INTERPOLATED LINE THROUGH
C XI(J1),ETA(J1) ...... XI(J),ETA(J) ...... XI(J1),ETA(J1)
C
   17	JJ = J1BIG
	IF (J1BIG.NE.1) JJ = J1BIG+1
CC	KNT = 0
	IP = 3
	DO 510 KK = JJ,J
CC		KNT = KNT+1
		X1=DX*(XI(KK)-XM)
		Y1=DY*(ETA(KK)-YM)
		CALL PLOT(X1,Y1,IP)
		IP = 2
CC		XX(KNT) = XI(KK)
CC		YY(KNT) = ETA(KK)
  510	CONTINUE
	X1=DX*(XI(JJ)-XM)
	Y1=DY*(ETA(JJ)-YM)
	CALL PLOT(X1,Y1,IP)
CC	XX(KNT+1) = XX(1)
CC	YY(KNT+1) = YY(1)
CC	NPOINT = KNT+1
CC	CALL CNTCRV (XX(1),YY(1),NPOINT,ZCON)
C
C (P)
C RESET J.  ESTABLISH THE NEXT CONTOUR LINE FOR REMAINING POINTS
C OR QUIT THE PROCEDURE IF NO MORE POINTS REMAIN.
C
	J = J1BIG - 1
	IF (J.GT.0) J=J+1
	IF (J) 800,800,11
  800	RETURN
	END
C
C
	SUBROUTINE TRIANGC (X,Y,NN,T,NZZ,M,P,B,NW,L,NZ,E,BE,TE)
C
C	COSMIC CONTOUR ROUTINES ARC-11441
C	MODIFIED  D.G. LONG, JPL  FEB 1986
C     ------------------------------------------------------------------
C
C A SET OF N DATA POINTS ARE KNOWN (X(I),Y(I),I=1,N)  THEY ARE TO
C BE CONNECTED BY LINES TO FORM A SET OF TRIANGLES. THE FINAL
C TRIANGULATION ESTABLISHES A CONVEX POLYGON DEFINED BY LINKED
C LISTS OF EDGE NUMBERS, END POINTS AND BOUNDARY EDGES.
C
C     INPUT
C       X   = ARRAY OF ABSCISSAS
C       Y   = ARRAY OF ORDINATES
C       NN  = NUMBER OF POINTS IN X AND Y
C		IF NN < 0 THEN BE AND TE ARRAYS ARE NOT USED
C	NW  = DIMENSION OF WORKING ARRAY
C	P,B = INTEGER WORKING ARRAYS P(NW),B(NW)
C	NZ  = DIMENSION OF INDEXING ARRAYS E,BE,TE
C	NZZ = DIMENSION OF INDEXING ARRAY T
C
C     OUTPUT
C       L   = NUMBER OF EDGES LISTED IN E, BE AND TE
C	T   = INDICES OF AJACENT TRIANGLE EDGES -- T(NZZ,3)
C  IF NN>0 THEN
C       E   = INTEGER LIST OF INDICES OF EACH TRIANGLE EDGE -- E(NZ,2)
C       BE  = 1 IF I OF E IS A BOUNDARY EDGE -- BE(NZ)
C       TE  = INDICES OF NEIGHBORING EDGES FOR EACH TRIANGLE -- TE(NZ,4)
C
C     LOCAL VARIABLES
C       P   = INDICES OF POINTS OUTSIDE THE BOUNDARY
C       J   = NO. OF VALUES IN LIST P
C       B   = INDEX OF POINTS ON THE BOUNDARY ..  INORDER
C       T   = INDICES OF ADJACENT TRIANGLE EDGES
C       K   = NO. OF POINTS LISTED IN ARRAY B
C       M   = NO. OF ROWS USED IN ARRAY T
C       X   = ARRAY OF SCALED X DATA
C       Y   = ARRAY OF SCALED Y DATA
C
C	EXTERNAL CALLS:
C		MIDC (INTEGER FUNCTION)
C
C     ------------------------------------------------------------------
C
	IMPLICIT INTEGER (P,B)
	INTEGER  T,TE,E
C
	DIMENSION  X(1),Y(1)
	DIMENSION  P(NW),B(NW)
	DIMENSION  T(NZZ,3)
C
	DIMENSION  E(NZ,2),BE(NZ),TE(NZ,4)
C
C DOUBLE PRECISION SPECIFICATION STATEMENTS
C
	DOUBLE PRECISION TERM,DCOMP,D,D1,S,TC
	DOUBLE PRECISION XP1,X21,YP1,Y21,XP2,X12,YP2,Y12,X1P,Y1P,X2P,Y2P
C
C (A)
C THE PROCEDURE BEGINS WITH NO BOUNDARY, NO EDGES, AND
C ALL X-Y DATA POINTS UNDER CONSIDERATION
C INITIALIZE LOCAL VARIABLES.
C
	N = IABS(NN)
	J = N
	K = 0
	L = 0
	M = 0
	KKNT = 0
	DO 100 JCNT=1,J
  100		P(JCNT) = JCNT
C
C (B)
C BEGIN BY TAKING THE LAST PAIR OF POINTS (X,Y(J)) IN THE
C LIST TO BE THE FIRST BOUNDARY POINT
C
	B(1) = J
	J = J-1
C
C (C)
C FROM THE REMAINING POINTS, FIND THE POINT NEAREST THE FIRST
C
	I2 = 1
	I1 = B(1)
	DMIN = (X(I1)-X(1))**2 + (Y(I1)-Y(1))**2
	DO 270 J1=2,J
		DST = (X(I1)-X(J1))**2 + (Y(I1)-Y(J1))**2
		IF (DST.GE.DMIN) GOTO 270
		I2 = J1
		DMIN = DST
  270	CONTINUE
C
C (D)
C NOW B(1) TO B(I2) IS THE FIRST EDGE.
C THERE IS ONE EDGE AND TWO BOUNDARY POINTS.
C
	J = J-1
	IF (I2.GT.J) GOTO 275
	DO 274 JCNT=I2,J
		P(JCNT) = P(JCNT+1)
  274	CONTINUE
  275	K = 2
	B(2) = I2
	L = 1
	IF (NN.GT.0) THEN
		E(1,1) = MIN0(B(1),B(2))
		E(1,2) = MAX0(B(1),B(2))
	ENDIF
C
C (E)
C NOW BEGIN CIRCLING AROUND THE BOUNDARY OF THE POLYGON,
C CONSIDERING, IN ORDER, EACH BOUNDARY EDGE.  MAINTAIN THE
C FOLLOWING INDICES -
C  K1 = B ARRAY INDEX OF THE CURRENT EDGE - POINT 1
C  K2 = B ARRAY INDEX OF THE CURRENT EDGE - POINT 2
C  B1,B2 = INDICES OF BOUNDARY POINT COORDINATES
C
	K1 = 0
	KT = 0
   11	K1 = K1+1
	IF (K1.GT.K) K1=1
   12	K2 = K1+1
	IF (K2.GT.K) K2=1
	B1 = B(K1)
	B2 = B(K2)
	KT = KT+1
C
C (F)
C CONSIDER THE BOUNDARY EDGE FROM B1 TO B2.  FOR ALL POINTS NOT
C YET TRIANGULATED (THE J POINTS REMAINING IN P), FIND THE
C POINT THAT, WHEN TRIANGULATED WITH B1,B2, MINIMIZES THE LENGTH
C OF THE TWO NEW EDGES TO BE DRAWN.
C
	D1 = 0.
	J1 = 0
	BFLAG = 0
	IF (J.EQ.0) GOTO 6
	DO 1 LJ=1,J
		PJ = P(LJ)
		TERM = (Y(PJ)-Y(B1))*(X(B2)-X(B1))-
     $			(X(PJ)-X(B1))*(Y(B2)-Y(B1))
		IF (TERM.LE.0.) GOTO 1
		D = SQRT((X(PJ)-X(B1))**2+(Y(PJ)-Y(B1))**2)
     $			+SQRT((X(PJ)-X(B2))**2+(Y(PJ)-Y(B2))**2)
		IF (J1.NE.0.AND.D1.LT.D) GOTO 1
		J1 = LJ
		D1 = D
    1	CONTINUE
C
C (G)
C IF LESS THAN THREE EDGES EXIST (NO TRIANGLE DEFINED YET),
C THEN THERE ARE NO ADJACENT BOUNDARY POINTS TO BE CONSIDERED.
C SO GO TO SECTION J.
C
	IF (K.LE.3) GOTO 3
C
C (H)
C CONSIDER THE ADJACENT BOUNDARY POINT OF THE NEXT EDGE OF THE
C POLYGON.  CALL ITS INDEX NUMBER K3 AND SEE IF ITS CLOSER TO
C THE CURRENT EDGE THAN P(J1).
C
    6	K3 = K2+1
	IF (K3.GT.K) K3=1
	PK3 = B(K3)
	TERM = (Y(PK3)-Y(B1))*(X(B2)-X(B1))-(X(PK3)-X(B1))*(Y(B2)-Y(B1))
	IF (TERM.LE.0.) GOTO 2
	D = SQRT((X(PK3)-X(B1))**2+(Y(PK3)-Y(B1))**2)
     $		+SQRT((X(PK3)-X(B2))**2+(Y(PK3)-Y(B2))**2)
	IF (J1.NE.0.AND.D1.LT.D) GOTO 2
	J1 = K3
	D1 = D
	BFLAG = 1
C
C (I)
C CONSIDER THE ADJACENT BOUNDARY POINT OF THE PREVIOUS EDGE OF
C THE POLYGON.  CALL ITS INDEX NUMBER K0 AND SEE IF ITS CLOSER
C TO THE CURRENT EDGE THAN P(J1) AND B(K3).
C
    2	CONTINUE
	K0 = K1-1
	IF (K0.LT.1) K0=K
	PK0 = B(K0)
	TERM = (Y(PK0)-Y(B1))*(X(B2)-X(B1))-(X(PK0)-X(B1))*(Y(B2)-Y(B1))
	IF (TERM.LE.0.) GOTO 3
	D = SQRT((X(PK0)-X(B1))**2+(Y(PK0)-Y(B1))**2)
     $		+SQRT((X(PK0)-X(B2))**2+(Y(PK0)-Y(B2))**2)
	IF (J1.NE.0.AND.D1.LT.D) GOTO 3
	J1 = K0
	D1 = D
	BFLAG = -1
    3	CONTINUE
C
C (J)
C SKIP THE NEXT SECTION IF J1 IS STILL ZERO, SINCE A CANDIDATE
C POINT FOR TRIANGULATION WITH EDGE B1,B2 WAS NOT FOUND.
C
	IF (J1.EQ.0) GOTO 9
C
C (K,L)
C IF THE SEARCH FOR A CANDIDATE POINT HAS ALREADY CONSIDERED EACH
C BOUNDARY EDGE AT LEAST ONCE (KT.GT.K) OR IF THE BOUNDARY IS
C BEING CHECKED FOR CONCAVE EDGES (J=0), THEN THE NEXT SECTION
C (SECTION M) CAN BE OMMITTED.
C
	IF (KT.GT.K.OR.J.EQ.0) GOTO 9
C
C (M)
C AT THIS POINT THE USER MAY INSERT ANY ADDITIONAL CONSTRAINT
C ON THE TRIANGLE TO BE FORMED BY THE POINT PJ1.  IF THE
C CANDIDATE TRIANGLE FAILS THE TEST, IT IS DELETED FROM
C CONSIDERATION BY SETTING THE VARIABLE J1 TO ZERO.
C
    9	CONTINUE
C
C (N,O)
C THE NEXT PROCEDURE CHECKS ALL BOUNDARY EDGES OF THE POLYGON
C FOR INTERSECTION WITH THE CANDIDATE TRIANGLE.  IF ANY EXISTING
C BOUNDARY EDGE INTERSECTS ANY OF THE EDGES TO BE FORMED BY THE
C CANDIDATE TRIANGLE, THEN THE CANDIDATE POINT IS REJECTED.  IF
C BFLAG IS NOT ZERO, THEN THE EDGE DEFINED BY J1=K0 OR J1=K3 IS
C EXEMPT FORM THIS TEST.
C
C IF THERE ARE THREE OR LESS EXISTING BOUNDARY EDGES OR IF
C J1 HAS BEEN SET TO ZERO, THIS TEST IS OMMITTED.
C
	IF (K.LE.3.OR.J1.EQ.0) GOTO 7
	IF (BFLAG.EQ.0) NQ = P(J1)
	IF (BFLAG.EQ.1) NQ = B(K3)
	IF (BFLAG.EQ.-1) NQ = B(K0)
	DO 108 KCNT=1,K
		IF (KCNT.EQ.K1) GOTO 108
		KN = KCNT+1
		IF (KCNT.EQ.K) KN=1
		IF (BFLAG.EQ.-1.AND.(KCNT.EQ.K0.OR.KN.EQ.K0)) GOTO 108
		IF (BFLAG.EQ. 1.AND.(KCNT.EQ.K3.OR.KN.EQ.K3)) GOTO 108
		P1 = B(KCNT)
		P2 = B(KN)
		DO 8 JCNT=1,2
			IF (JCNT.EQ.1.AND.(BFLAG.EQ.0.OR.BFLAG.EQ.1)
     $				.AND.KCNT.EQ.K0) GOTO 108
			IF (JCNT.EQ.2.AND.(BFLAG.EQ.0.OR.BFLAG.EQ.-1)
     $				.AND.KCNT.EQ.K2) GOTO 108
			BJ = B1
			IF (JCNT.EQ.2) BJ=B2
			XQB = X(NQ)-X(BJ)
			YQB = Y(NQ)-Y(BJ)
			X12 = X(P1)-X(P2)
			Y12 = Y(P1)-Y(P2)
			D = XQB*Y12-YQB*X12
			IF (D.EQ.0.) GOTO 8
			X1B = X(P1)-X(BJ)
			Y1B = Y(P1)-Y(BJ)
			S = (X1B*Y12-Y1B*X12)/D
			IF (S.LT.0..OR.S.GT.1.) GOTO 8
			TC = (XQB*Y1B-YQB*X1B)/D
			IF (TC.LT.0.OR.TC.GT.1.) GOTO 8
			J1 = 0
			GOTO 7
    8		CONTINUE
  108 	     CONTINUE
    7	CONTINUE
C
C (P,Q)
C IF J1 IS ZERO, THEN THE CANDIDATE POINT DID NOT PASS THE ABOVE
C TESTS OR NO POINT WAS FOUND.  IF BFLAG IS NOT ZERO, THEN A
C POINT ON THE BOUNDARY WAS FOUND.
C
	IF (J1.EQ.0) GOTO 10
	IF (BFLAG) 150,160,4
C
C THE TRIANGULATED POINT IS OUTSIDE THE BOUNDARY.  ESTABLISH TWO
C NEW EDGES, A NEW BOUNDARY POINT AND DELETE ONE POINT FROM
C OUTSIDE THE BOUNDARY.
C
  160	IF (NN.GT.0) THEN
		E(L+1,1) = MIN0(P(J1),B(K1))
		E(L+1,2) = MAX0(P(J1),B(K1))
		E(L+2,1) = MIN0(P(J1),B(K2))
		E(L+2,2) = MAX0(P(J1),B(K2))
	ENDIF
	KT = 0
	L = L+2
	M = M+1
	T(M,1) = MIN0(P(J1),B(K1),B(K2))
	T(M,2) = MIDC(P(J1),B(K1),B(K2))
	T(M,3) = MAX0(P(J1),B(K1),B(K2))
	IF (K1.EQ.K) GOTO 140
	 KM = K
	 KP1 = K1+1
  147	 B(KM+1) = B(KM)
	 KM = KM-1
	 IF (KM.GE.KP1) GOTO 147
  140	B(K1+1) = P(J1)
	K = K+1
	J = J-1
	IF (J1.GT.J) GOTO 10
	DO 144 JCNT=J1,J
  144		P(JCNT) = P(JCNT+1)
	GOTO 10
C
C (S)
C THE TRIANGULATED POINT IS THE NEXT POINT ON THE BOUNDARY.
C ESTABLISH ONE NEW EDGE (FROM B(K1) TO B(K3)), ONE NEW
C TRIANGLE (FROM B(K1) TO B(K2) TO B(K3)), AND DELETE ONE POINT
C FROM THE BOUNDARY (B(K2)).
C
    4	IF (NN.GT.0) THEN
		E(L+1,1) = MIN0(B(K3),B(K1))
		E(L+1,2) = MAX0(B(K3),B(K1))
	ENDIF
	KK = 0
	KKNT = 0
	KT = 0
	L=L+1
	K=K-1
	M = M+1
	T(M,1) = MIN0(B(K1),B(K2),B(K3))
	T(M,2) = MIDC(B(K1),B(K2),B(K3))
	T(M,3) = MAX0(B(K1),B(K2),B(K3))
	IF (K2.GT.K) GOTO 155
	DO 151 KCNT=K2,K
  151		B(KCNT) = B(KCNT+1)
  155	IF (K2.EQ.1) K1=K1-1
	GOTO 10
C
C (R)
C THE TRIANGULATED POINT IS THE PREVIOUS POINT ON THE BOUNDARY.
C ESTABLISH A NEW EDGE (FROM B(K0) TO B(K2)), ONE NEW TRIANGLE
C (FROM B(K0) TO B(K1) TO B(K2)), AND DELETE ONE POINT FROM THE
C BOUNDARY (B(K1))
C
  150	IF (NN.GT.0) THEN
		E(L+1,1) = MIN0(B(K0),B(K2))
		E(L+1,2) = MAX0(B(K0),B(K2))
	ENDIF
	KK = 0
	KKNT = 0
	KT = 0
	L = L+1
	K = K-1
	M = M+1
	T(M,1) = MIN0(B(K0),B(K1),B(K2))
	T(M,2) = MIDC(B(K0),B(K1),B(K2))
	T(M,3) = MAX0(B(K0),B(K1),B(K2))
	IF (K1.GT.K) GOTO 157
	DO 158 KCNT=K1,K
  158		B(KCNT) = B(KCNT+1)
  157	K1 = K1-1
	IF (K1.LT.1) K1=K
C
C (T)
C IF THERE ARE ANY POINTS REMAINING OUTSIDE THE BOUNDARY, THEN
C REPEAT THE PROCEDURE FOR THE NEXT EDGE.
C
   10	IF (J.GT.0.AND.J1.NE.0) GOTO 12
	IF (J.GT.0) GOTO 11
C
C (U,V,W)
C ALL POINTS HAVE BEEN TRIANGULATED.  CHECK THAT ALL BOUNDARY
C EDGES FORM A CONCAVE POLYGON.
C
	IF (KK.NE.0) GOTO 55
	KK = 1
	KL = 0
   55	KKNT = KKNT+1
	IF (KKNT.GT.N) GOTO 170
    5	KL = KL+1
	K2 = KL+1
	IF (K2.GT.K) K2=1
	K1 = KL-1
	IF (K1.LT.1) K1=K
	PKL = B(KL)
	B1 = B(K1)
	B2 = B(K2)
	TERM = (Y(PKL)-Y(B1))*(X(B2)-X(B1))-(X(PKL)-X(B1))*(Y(B2)-Y(B1))
	IF (TERM.LT.0.) GOTO 11
	IF (KL.LT.K) GOTO 5
C
C (X)
C THE TRIANGULATION IS COMPLETE AND HAS BEEN CHECKED FOR A
C CONCAVE BOUNDARY.  NOW IDENTIFY THE BOUNDARY EDGES.
C
	IF (NN.LT.0) RETURN
  170	DO 23 LCNT=1,L
	BE(LCNT) = 0
	KL = 0
   21	KL = KL+1
	IF (E(LCNT,1).NE.B(KL)) GOTO 22	!%
	K1 = KL+1
	IF (K1.GT.K) K1=1
	IF (E(LCNT,2).NE.B(K1)) GOTO 162!%
	BE(LCNT) = 1
	GOTO 23
  162	K1 = KL-1
	IF (K1.LT.1) K1=K
	IF (E(LCNT,2).NE.B(K1)) GOTO 22 !%
	IF (NN.GT.0) BE(LCNT) = 1
	GOTO 23
   22	IF (KL.LT.K) GOTO 21
   23	CONTINUE
C
C (Y)
C FINALLY, ESTABLISH THE INDICES OF ADJACENT EDGES FOR EACH
C EDGE IN THE TRIANGULATION.  EACH BOUNDARY EDGE WILL HAVE TWO
C ADJACENT EDGES - EACH INTERIOR EDGE WILL HAVE FOUR.
C
	DO 190 LL  =1,4
		DO 190 LCNT=1,L
  190			TE(LCNT,LL) = 0
	DO 191 MCNT=1,M
		DO 192 LL=1,L
		  IF (E(LL,1).EQ.T(MCNT,1).AND.E(LL,2).EQ.T(MCNT,2)) L1=LL
		  IF (E(LL,1).EQ.T(MCNT,2).AND.E(LL,2).EQ.T(MCNT,3)) L2=LL
		  IF (E(LL,1).EQ.T(MCNT,1).AND.E(LL,2).EQ.T(MCNT,3)) L3=LL
  192		CONTINUE
		LAMBDA = 0
		IF (TE(L1,1).NE.0) LAMBDA=2
		TE(L1,LAMBDA+1) = L2
		TE(L1,LAMBDA+2) = L3
		LAMBDA = 0
		IF (TE(L2,1).NE.0) LAMBDA=2
		TE(L2,LAMBDA+1) = L1
		TE(L2,LAMBDA+2) = L3
		LAMBDA = 0
		IF (TE(L3,1).NE.0) LAMBDA = 2
		TE(L3,LAMBDA+1) = L1
		TE(L3,LAMBDA+2) = L2
  191	CONTINUE
C
	RETURN
	END
C
	FUNCTION MIDC(I,J,K)
C
C THIS FUNCTION SUBPROGRAM IS USED BY THE TRIANGULATION ALGORITHM
C TO FIND THE MIDDLE VALUE OF THE THREE INTEGER ARGUMENTS (THE
C VALUE WHICH IS NEITHER A MINIMUM OR A MAXIMUM).  I, J AND K ARE
C ARE ASSUMED TO BE DISCRETE VALUES WITH NO TWO EQUAL.
C
	IF (J.LT.I.AND.I.LT.K) GOTO 100
	IF (K.LT.I.AND.I.LT.J) GOTO 100
	IF (I.LT.J.AND.J.LT.K) GOTO 200
	IF (K.LT.J.AND.J.LT.I) GOTO 200
	MIDC = K
	RETURN
  100	MIDC = I
	RETURN
  200	MIDC = J
	RETURN
	END
C
C **************************************************************************
C
	SUBROUTINE TRIG3DH(X,Y,Z,N,PSI,PHI,THETA,
     $		DV,XH,YH,ZH,IFLAG,IAXIS,
     $		XT,NXT,TX,TSX,PGX,FSX,
     $		YT,NYT,TY,TSY,PGY,FSY,
     $		ZT,NZT,TZ,TSZ,PGZ,FSZ,
     $		XAST,XAEN,YAST,YAEN,ZAST,ZAEN,ICOL)
C
C	CREATED BY D. LONG     FEB, 1986	AT JPL
C
C	ROUTINE TO PLOT DATA SPECIFIED AS (X,Y,Z) TRIPLES AS 3-D SURFACE
C	BY TRIANGULATING AND PLOTTING WITH HIDDEN LINE REMOVAL
C	USING THE INIT3DH ROUTINES.
C
C	X,Y,Z	 (R) ARRAYS CONTAINING DATA POINTS
C	N	 (I) NUMBER OF DATA POINTS
C	PSI,PHI,THETA (R) ANGLES (IN DEGREES) OF INIT3DH
C	DV	 (R) PERSPECTIVE FACTOR OF INIT3DH (9999=NO PERSPECTIVE)
C		    DV < 0 DO NOT INITIALIZE INIT3DH
C	XH,YH,ZH (R) LENGTH OF EACH AXIS (INCHES)
C	IFLAG	(I) PLOT FLAG
C		>1000	REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0	REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT)
C		= 0	CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0	SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C 	(ONE'S DIGIT) = 2 USE RAM TEK COLOR CONTROL ARRAY
C		      = 1 DO NOT USE RAM TEK COLOR ARRAY
C	(TEN'S DIGIT) = 0 NO HIDDEN LINE REMOVAL FOR AXIS
C		      = 1 HIDDEN LINE REMOVAL OF AXIS
C		      = 2 NO HIDDEN LINES AT ALL (SKETCH USED)
C		      = 3 NO HIDDEN LINES AT ALL (SKETCH NOT USED)
C 	(HUNDRED'S)   = 0 ASK WHICH SCREEN DEVICE TO USE
C		      <>0 SCREEN DEVICE CODE NUMBER
C 	IAXIS	(I) AXIS OPTION FLAG
C			< 0 DO NOT PLOT AXES
C 	(ONE'S DIGIT)	= 1 PLOT Z AXIS USING COMPUTED SCALE
C			= 2 PLOT Z AXIS USING COMPUTED, SMOOTHED SCALE
C			= 3 PLOT Z AXIS USING INPUT SCALE--VARIABLES REQUIRED
C			= 4 PLOT Z AXIS USING INPUT, SMOOTHED SCALE--VARIABLES REQUIRED
C 	(TEN'S DIGIT)	= 1 PLOT Y AXIS USING COMPUTED SCALE
C			= 2 PLOT Y AXIS USING COMPUTED, SMOOTHED SCALE
C			= 3 PLOT Y AXIS USING INPUT SCALE--VARIABLES REQUIRED
C			= 4 PLOT Y AXIS USING INPUT, SMOOTHED SCALE--VARIABLES REQUIRED
C 	(HUNDRED'S)	= 1 PLOT X AXIS USING COMPUTED SCALE
C			= 2 PLOT X AXIS USING COMPUTED, SMOOTHED SCALE
C			= 3 PLOT X AXIS USING INPUT SCALE--VARIABLES REQUIRED
C			= 4 PLOT X AXIS USING INPUT, SMOOTHED SCALE--VARIABLES REQUIRED
C 	XT,YT,ZT	(B) STRINGS FOR AXIS TITLES
C 	NXT,NYT,NZT	(I) LENGTH OF AXIS TITLES
C			     IF ZERO THEN AXIS NOT PLOTTED
C FOLLOWING ONLY REQUIRED IF IAXIS > 0
C	TX,TY,TZ	(R) NUMBER OF TICK MARKS (SEE AXIS3DH)
C	TSX,TSY,TSZ	(R) SIZE OF TITLE AND NUMBERS OF AXIS
C			  IF LESS THAN ZERO DO NOT AUTO-SCALE BY (x10^POWER)
C	PGX,PGY,PGZ	(R) ROTATION OF TEXT AND TICKS AROUND AXIS (DEG)
C	FSX,FSY,FSZ	(R) NUMBER FORMAT OF AXISES (SEE AXIS3DH)
C 	XAST,YAST,ZAST	(R) AXIS START VALUES
C 	XAEN,YAEN,ZAEN	(R) AXIS END VALUES
C	ICOL		(I) COLOR INDEX TABLE (REQUIRED IF 1'S DIGIT OF IFLAG=2
C				ICOL(1) AXIS COLOR
C				ICOL(2) SURFACE COLOR
C			    NOTE: WHEN BOTH AXIS AND SURFACE USE SKETCH
C			    BOTH ARE PLOTTED USING SURFACE COLOR
C
	DIMENSION X(1),Y(1),Z(1),PAS(2),ICOL(2)
	LOGICAL REPEAT/.FALSE./,HIDDEN
	INTEGER XT(1),YT(1),ZT(1)
C
	DIMENSION XX(4),YY(4),ZZ(4)
	DIMENSION ITT(2000,3),IP(1000),IB(1000)
	DATA NZZ/2000/,NW/1000/
C
	PARAMETER (IWRKSZE=100000)
	COMMON /GO/ISIZE,WORK(IWRKSZE)
C
	IF (ISIZE.LT.IWRKSZE) ISIZE=IWRKSZE	! INCREASE WORKSPACE IF NEEDED
C
	IF (N.GT.NW) THEN
		CALL CTERM(1)
		WRITE (*,3010)
3010		FORMAT(' *** WARNING: INSUFFICIENT WORKING SPACE')
		CALL CTERM(-1)
	ENDIF
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	XLEN=ABS(XH)
	YLEN=ABS(YH)
	ZLEN=ABS(ZH)
	IZAXIS=MOD(IABS(IAXIS),10)
	IYAXIS=MOD(IABS(IAXIS),100)/10
	IXAXIS=MOD(IABS(IAXIS),1000)/100
	IF (IXAXIS.GT.2) THEN
		XMAX=XAEN
		XMIN=XAST
	ELSE
		XMAX=X(1)
		XMIN=X(1)
		DO 16 I=1,N
			XMAX=AMAX1(XMAX,X(I))
			XMIN=AMIN1(XMIN,X(I))
16		CONTINUE
	ENDIF
	IF (IXAXIS.EQ.2.OR.IXAXIS.EQ.4) THEN
		PAS(1)=XMAX
		PAS(2)=XMIN
		CALL SCALE(PAS,XLEN,2,1,1,XMIN,DAA)
		XMAX=YLEN*DAA+XMIN
	ENDIF
C
	IF (IYAXIS.GT.2) THEN
		YMAX=YAEN
		YMIN=YAST
	ELSE
		YMAX=Y(1)
		YMIN=Y(1)
		DO 17 I=1,N
			YMAX=AMAX1(YMAX,Y(I))
			YMIN=AMIN1(YMIN,Y(I))
17		CONTINUE
	ENDIF
	IF (IYAXIS.EQ.2.OR.IYAXIS.EQ.4) THEN
		PAS(1)=YMAX
		PAS(2)=YMIN
		CALL SCALE(PAS,YLEN,2,1,1,YMIN,DAA)
		YMAX=YLEN*DAA+YMIN
	ENDIF
C
	IF (IZAXIS.GT.2) THEN
		ZMAX=ZAEN
		ZMIN=ZAST
	ELSE
		ZMAX=Z(1)
		ZMIN=Z(1)
		DO 18 I=1,N
			ZMAX=AMAX1(ZMAX,Z(I))
			ZMIN=AMIN1(ZMIN,Z(I))
18		CONTINUE
	ENDIF
	IF (IZAXIS.EQ.2.OR.IZAXIS.EQ.4) THEN
		PAS(1)=ZMAX
		PAS(2)=ZMIN
		CALL SCALE(PAS,ZLEN,2,1,1,ZMIN,DAA)
		ZMAX=ZLEN*DAA+ZMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=MOD(JF,10000)/100
		CALL FRAME(3,ILU,1.5,0.65,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL IN PLOT MODE
	ENDIF
C
	HIDDEN=.TRUE.
	IOPT=MOD(JF,100)/10
	IF (IOPT.EQ.3) THEN
		IOPT=0
		HIDDEN=.FALSE.
	ENDIF
	JF=MOD(JF,10)
C
C	INITIALIZE 3DH PACKAGE
C
      	MNE=3			! MAX NUMBER OF EDGES
      	MNH=0			! NO HOLES
	ICORE=(25+5*MNE+4*MNH)*NX*NY
	IF (IAXIS.NE.0) ICORE=ICORE+10000	! ESTIMATED AXISES SPACE
	IF (ICORE.GT.ISIZE) THEN
		CALL CTERM(1)
		WRITE (*,3011)
3011		FORMAT('*** WARNING: COMMON BLOCK /GO/ MAY BE TOO SMALL')
		CALL CTERM(-1)
	ENDIF
	SF=1.0			! SCALE FACTOR
	IF (IOPT.EQ.2) SF=-1.0	! NO HIDDEN LINE REMOVAL
	IF (DV.GT.0.) CALL INIT3DH(0.,0.,0.,PSI,PHI,THETA,SF,DV,MNE,MNH)
C
	IF (IAXIS.NE.0) THEN		!PLOT AXIS LABELS
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(1)),0.,0) !RAM TEK COLOR
		IF (NXT.GT.0) THEN 	!PLOT X AXIS
			CALL AXIS3DH(0.,0.,0.,0.,0.,PGX,XT,
     $				NXT,XLEN,XMIN,XMAX,
     $				TX,TSX,FSX,IOPT)
		ENDIF
		IF (NYT.GT.0) THEN 	!PLOT Y AXIS
			CALL AXIS3DH(0.,0.,0.,0.,90.,PGY,YT,
     $				NYT,YLEN,YMIN,YMAX,
     $				TY,TSY,FSY,IOPT)
		ENDIF
		IF (NZT.GT.0) THEN 	!PLOT Z AXIS
			CALL AXIS3DH(0.,0.,0.,90.,0.,PGZ,ZT,
     $				NZT,-ZLEN,ZMIN,ZMAX,
     $				TZ,TSZ,FSZ,IOPT)
		ENDIF
	ENDIF
	IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(2)),0.,0) !RAM TEK COLOR
C
C	BEGIN SURFACE PLOTTING
C
	DX=1.0
	DY=1.0
	DZ=1.0
	IF (XMAX-XMIN.NE.0.0) DX=XLEN/(XMAX-XMIN)
	IF (YMAX-YMIN.NE.0.0) DY=YLEN/(YMAX-YMIN)
	IF (ZMAX-ZMIN.NE.0.0) DZ=ZLEN/(ZMAX-ZMIN)
C
	CALL TRIANGC(X,Y,-N,ITT,NZZ,M,IP,IB,NW,L,NZZ,IE,IBE,ITE)
	DO 50 I=1,M
		DO 15 J=1,3
			K1=ITT(I,J)
			XX(J)=(X(K1)-XMIN)*DX
			YY(J)=(Y(K1)-YMIN)*DY
			ZZ(J)=(Z(K1)-ZMIN)*DY
15		CONTINUE
		XX(4)=XX(1)
		YY(4)=YY(1)
		ZZ(4)=ZZ(1)
		NC=4
		IERR=0
		IF (I.EQ.M) IERR=1
		IF (HIDDEN) THEN
			CALL SKETCH(XX,YY,ZZ,NC,IERR)
		ELSE
			IERR=3
			DO 25 JF=1,NC
			  CALL PLT3DH(XX(JF),YY(JF),ZZ(JF),IERR)
			  IERR=2
25			CONTINUE
		ENDIF
		IF (IERR.LT.0) THEN
			CALL CTERM(1)
			WRITE(*,3012) IERR
3012			FORMAT(' *** WARNING: SKETCH ERROR HAS OCCURED',I5)
			CALL CTERM(-1)
		ENDIF
 50	CONTINUE
	IF (.NOT.HIDDEN) CALL PLT3DH(XX(JF),YY(JF),ZZ(JF),3)	!PEN UP
C
C	FINISH ROUTINE
C
999	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL IN TEXT MODE
	ENDIF
	RETURN
	END
C
C **************************************************************************
C
	SUBROUTINE HIST3D(A,INX,INY,NX,NY,PSI,PHI,THETA,
     $		DV,ZMIN,XH,YH,ZH,IFLAG,IAXIS,
     $		XT,NXT,XASTART,XAEND,TX,TSX,PGX,FSX,
     $		YT,NYT,YASTART,YAEND,TY,TSY,PGY,FSY,
     $		ZT,NZT,TZ,TSZ,PGZ,FSZ,AMININ,AMAXIN,ICOL)
C
C	CREATED BY D. LONG     FEB, 1986	AT JPL
C
C	ROUTINE TO PLOT DATA IN 3-D COLUMN FORM WITH HIDDEN LINE REMOVAL
C	USING THE INIT3DH ROUTINES.
C
C	A	(R) ARRAY A(INX,INY) CONTAINING Z AXIS DATA
C	INX,INY (I) DIMENSION OF A ARRAY
C	NX,NY	(I) INDICATING SIZE OF A ARRAY TO PLOT
C	PSI,PHI,THETA (R) ANGLES (IN DEGREES) OF INIT3DH
C	DV	(R) PERSPECTIVE FACTOR OF INIT3DH (9999 = NO PERSPECTIVE)
C		    DV < 0 DO NOT INITIALIZE INIT3DH
C	ZMIN	(R) REFERENCE PLANE VALUE FOR Z AXIS
C	XH,YH,ZH (R) LENGTH OF EACH AXIS (INCHES)
C	IFLAG	(I) PLOT FLAG
C		>1000	REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0	REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT
C		= 0	CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0	SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C 	(ONE'S DIGIT) = 2 USE RAM TEK COLOR CONTROL ARRAY
C		      = 1 DO NOT USE RAM TEK COLOR ARRAY
C	(TEN'S DIGIT) = 0 NO HIDDEN LINE REMOVAL FOR AXIS
C		      = 1 HIDDEN LINE REMOVAL OF AXIS
C		      = 2 NO HIDDEN LINES (SKETCH STILL USED)
C		      = 3 NO HIDDEN LINES (SKETCH NOT USED)
C 	(HUNDRED'S)   = 0 ASK WHICH SCREEN DEVICE TO USE
C		      <>0 SCREEN DEVICE CODE NUMBER
C 	IAXIS	(I) AXIS OPTION FLAG
C			= 0 DO NOT PLOT AXIS--FOLLOWING VARIABLES NOT REQUIRED
C			< 0 PLOT AXIS, USE INPUT Z AXIS SCALE
C			> 0 PLOT AXIS, USE COMPUTED Z AXIS SCALE
C 	(ONE'S DIGIT)	= 1 PLOT AXIS, Y AXIS SCALE--VARIABLES REQUIRED
C			= 2 PLOT AXIS, AUTO SCALE Y AXIS--VARIABLES REQUIRED
C 	XT,YT,ZT	(B) STRINGS FOR AXIS TITLES
C 	NXT,NYT,NZT	(I) LENGTH OF AXIS TITLES
C			     IF ZERO THEN THAT AXIS NOT PLOTTED
C 	XASTART,YASTART	(R) AXIS START VALUES
C 	XAEND,YAEND	(R) AXIS END VALUES
C FOLLOWING ONLY REQUIRED IF IAXIS<>0
C	TX,TY,TZ	(R) NUMBER OF TICK MARKS (SEE AXIS3DH)
C	TSX,TSY,TSZ	(R) SIZE OF TITLE AND NUMBERS OF AXIS
C			  IF LESS THAN ZERO DO NOT AUTO-SCALE BY (x10^POWER)
C	PGX,PGY,PGZ	(R) ROTATION OF TEXT AND TICKS AROUND AXIS (DEG)
C	FSX,FSY,FSZ	(R) NUMBER FORMAT OF AXISES (SEE AXIS3DH)
C 	AMININ,AMAXIN 	(R) ZAXIS SCALING FACTORS (ONLY NEEDED IF IAXIS < 0)
C	ICOL		(I) COLOR INDEX TABLE (REQUIRED IF 1'S DIGIT OF IFLAG=2)
C				ICOL(1) : AXIS COLOR
C				ICOL(2) : COLUMNS
C			    NOTE: IF SKETCH IS USED FOR BOTH AXIS AND COLUMNS
C			    THEN ENTIRE PLOT IS IN COLUMN COLOR
C
	DIMENSION A(INX,INY),PAS(2),ICOL(2)
	LOGICAL REPEAT/.FALSE./,HIDDEN
	INTEGER XT(1),YT(1),ZT(1)
C
	DIMENSION XX(5),YY(5),ZZ(5)
C
	PARAMETER (IWRKSZE=100000)
	COMMON /GO/ISIZE,WORK(IWRKSZE)
C
C	INSURE THAT WE HAVE AT LEAST A MINIMUM SIZE WORKSPACE
C
	IF (ISIZE.LT.IWRKSZE) ISIZE=IWRKSZE	! INCREASE WORKSPACE IF NEEDED
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
C	DETERMINE MAX/MIN Z VALUES
C
	IF (IAXIS.LT.0) THEN
		AMAX=AMAXIN
		AMIN=AMININ
	ELSE
		AMAX=ZMIN
		AMIN=ZMIN
		DO 16 IY=1,NY			!DETERMINE MAX,MIN ARRAY VALUES
			DO 16 IX=1,NX
				AMAX=AMAX1(AMAX,A(IX,IY))
				AMIN=AMIN1(AMIN,A(IX,IY))
16		CONTINUE
	ENDIF
C
C	SMOOTH Z VALUES IF DESIRED
C
	XLEN=ABS(XH)
	YLEN=ABS(YH)
	ZLEN=ABS(ZH)
	IF (MOD(IABS(IAXIS),10).EQ.2) THEN
		PAS(1)=AMAX
		PAS(2)=AMIN
		CALL SCALE(PAS,ZLEN,2,1,1,AMIN,DAA)
		AMAX=ZLEN*DAA+AMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=MOD(JF,10000)/100
		CALL FRAME(3,ILU,1.5,0.65,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL IN PLOT MODE
	ENDIF
C
	HIDDEN=.TRUE.
	IOPT=MOD(JF,10000)/10
	IF (IOPT.EQ.3) THEN
		IOPT=0
		HIDDEN=.FALSE.
	ENDIF
	JF=MOD(JF,10)
C
C	INITIALIZE 3DH PACKAGE
C
      	MNE=4			! MAX NUMBER OF EDGES
      	MNH=0			! NO HOLES
	ICORE=(25+5*MNE+4*MNH)*NX*NY*6
	IF (IAXIS.NE.0) ICORE=ICORE+10000	! ESTIMATED AXISES SPACE
	IF (ICORE.GT.ISIZE) THEN
		CALL CTERM(1)
		WRITE(*,3030)
3030		FORMAT(' *** WARNING: COMMON BLOCK /GO/ MAY BE TOO SMALL')
		CALL CTERM(-1)
	ENDIF
	SF=1.0			! SCALE FACTOR
	IF (IOPT.EQ.2) SF=-1.0	! NO HIDDEN LINE REMOVAL FOR SKETCH
	IF (DV.GT.0.) CALL INIT3DH(0.,0.,0.,PSI,PHI,THETA,SF,DV,MNE,MNH)
C
C	PLOT AXES IF DESIRED
C
	IF (IAXIS.NE.0) THEN		!PLOT AXIS LABELS
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(1)),0.,0) !RAM TEK COLOR
		IF (NXT.GT.0) THEN 	!PLOT X AXIS
			CALL AXIS3DH(0.,0.,0.,0.,0.,PGX,XT,
     $				NXT,XLEN,XASTART,XAEND,
     $				TX,TSX,FSX,IOPT)
		ENDIF
		IF (NYT.GT.0) THEN 	!PLOT Y AXIS
			CALL AXIS3DH(0.,0.,0.,0.,90.,PGY,YT,
     $				NYT,YLEN,YASTART,YAEND,
     $				TY,TSY,FSY,IOPT)
		ENDIF
		IF (NZT.GT.0) THEN 	!PLOT Z AXIS
			CALL AXIS3DH(0.,0.,0.,90.,0.,PGZ,ZT,
     $				NZT,-ZLEN,AMIN,AMAX,
     $				TZ,TSZ,FSZ,IOPT)
		ENDIF
	ENDIF
	IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(2)),0.,0) !RAM TEK COLOR
C
C	BEGIN SURFACE PLOTTING
C
	XXX=NX
	DX=XLEN/IABS(NX+1)
	DY=YLEN/IABS(NY+1)
	DZ=1.0
	IF (AMAX-AMIN.NE.0.0) DZ=ZLEN/(AMAX-AMIN)
	NC=5
	Z1=ZMIN
	IERR=0
C
	DO 50 IY=1,IABS(NY)		! FOR EACH Y CELL
		Y1=(IY-1)*DY
		Y2=IY*DY
		DO 50 IX=1,IABS(NX)	! FOR EACH X CELL
			X1=(IX-1)*DX
			X2=IX*DX
			Z2=(A(IX,IY)-AMIN)*DZ
			IERR=0
			IF (IX.EQ.IABS(NX).AND.IY.EQ.IABS(NY)) IERR=1
			IF(Z1.NE.Z2)CALL CUBE(X1,X2,Y1,Y2,Z1,Z2,IERR,HIDDEN)
			IF (IERR.LT.0) THEN
			 CALL CTERM(1)
			 WRITE(*,3031) IERR
3031		FORMAT(' *** WARNING: SKETCH ERROR HAS OCCURED',I3)
			 CALL CTERM(-1)
			ENDIF
 50	CONTINUE
C
C	FINISH ROUTINE
C
	IF (.NOT.HIDDEN) CALL PLT3DH(X1,Y1,Z1,3)	!PEN UP
999	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL TEXT MODE
	ENDIF
	RETURN
	END
C
	SUBROUTINE CUBE(X1,X2,Y1,Y2,Z1,Z2,IERR,HIDE)
	LOGICAL HIDE
C
C	PLOTS A 3D SOLID 'CUBE' USING SKETCH OR PLOT3DH
C
C	X1,Y1,Z1 (R) ONE CORNER OF CUBE
C	X2,Y2,Z2 (R) OPOSITE CORNER OF CUBE
C	IERR	 (I) OPTION/ERR FLAG FOR LAST CALL TO SKETCH
C	HIDE	 (L) IF TRUE USE SKETCH ELSE USE PLT3DH
C
	DIMENSION XX(5),YY(5),ZZ(5)
C
	IER=0
	NC=5
C
	XX(1)=X1
	XX(2)=X2
	XX(3)=X2
	XX(4)=X1
	XX(5)=X1
	YY(1)=Y1
	YY(2)=Y1
	YY(3)=Y2
	YY(4)=Y2
	YY(5)=Y1
	ZZ(1)=Z1
	ZZ(2)=Z1
	ZZ(3)=Z1
	ZZ(4)=Z1
	ZZ(5)=Z1
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IER)	! BOTTOM
	ELSE
		IP=3
		DO 1 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
1		CONTINUE
	ENDIF
	IF (IER.LT.0) GOTO 10
	XX(1)=X1
	XX(2)=X2
	XX(3)=X2
	XX(4)=X1
	XX(5)=X1
	YY(1)=Y1
	YY(2)=Y1
	YY(3)=Y2
	YY(4)=Y2
	YY(5)=Y1
	ZZ(1)=Z2
	ZZ(2)=Z2
	ZZ(3)=Z2
	ZZ(4)=Z2
	ZZ(5)=Z2
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IER)	! TOP
	ELSE
		IP=3
		DO 2 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
2		CONTINUE
	ENDIF
	IF (IER.LT.0) GOTO 10
	XX(1)=X1
	XX(2)=X1
	XX(3)=X1
	XX(4)=X1
	XX(5)=X1
	YY(1)=Y1
	YY(2)=Y1
	YY(3)=Y2
	YY(4)=Y2
	YY(5)=Y1
	ZZ(1)=Z1
	ZZ(2)=Z2
	ZZ(3)=Z2
	ZZ(4)=Z1
	ZZ(5)=Z1
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IER)	! LEFT SIDE
	ELSE
		IP=3
		DO 3 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
3		CONTINUE
	ENDIF
	IF (IER.LT.0) GOTO 10
	XX(1)=X2
	XX(2)=X2
	XX(3)=X2
	XX(4)=X2
	XX(5)=X2
	YY(1)=Y1
	YY(2)=Y1
	YY(3)=Y2
	YY(4)=Y2
	YY(5)=Y1
	ZZ(1)=Z1
	ZZ(2)=Z2
	ZZ(3)=Z2
	ZZ(4)=Z1
	ZZ(5)=Z1
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IER)	! RIGHT SIDE
	ELSE
		IP=3
		DO 4 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
4		CONTINUE
	ENDIF
	IF (IER.LT.0) GOTO 10
	XX(1)=X1
	XX(2)=X2
	XX(3)=X2
	XX(4)=X1
	XX(5)=X1
	YY(1)=Y1
	YY(2)=Y1
	YY(3)=Y1
	YY(4)=Y1
	YY(5)=Y1
	ZZ(1)=Z1	
	ZZ(2)=Z1
	ZZ(3)=Z2
	ZZ(4)=Z2
	ZZ(5)=Z1
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IER)	! FRONT
	ELSE
		IP=3
		DO 5 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
5		CONTINUE
	ENDIF
	IF (IER.LT.0) GOTO 10
	XX(1)=X1
	XX(2)=X2
	XX(3)=X2
	XX(4)=X1
	XX(5)=X1
	YY(1)=Y2
	YY(2)=Y2
	YY(3)=Y2
	YY(4)=Y2
	YY(5)=Y2
	ZZ(1)=Z1
	ZZ(2)=Z1
	ZZ(3)=Z2
	ZZ(4)=Z2
	ZZ(5)=Z1
	IF (HIDE) THEN
		CALL SKETCH(XX,YY,ZZ,NC,IERR)	! BACK
	ELSE
		IP=3
		DO 6 I=1,5
			CALL PLT3DH(XX(I),YY(I),ZZ(I),IP)
			IP=2
6		CONTINUE
	ENDIF
10	IERR=IER
	RETURN
	END
C
C **************************************************************************
C
	SUBROUTINE T3DH(A,INX,INY,NX,NY,PSI,PHI,THETA,
     $		DV,XH,YH,ZH,IFLAG,IAXIS,
     $		XT,NXT,XASTART,XAEND,TX,TSX,PGX,FSX,
     $		YT,NYT,YASTART,YAEND,TY,TSY,PGY,FSY,
     $		ZT,NZT,TZ,TSZ,PGZ,FSZ,AMININ,AMAXIN,ICOL)
C
C	CREATED BY D. LONG     FEB, 1986	AT JPL
C
C	ROUTINE TO PLOT DATA IN A TRIANGULAR 3-D MESH FORM WITH HIDDEN LINE
C	REMOVAL USING THE INIT3DH ROUTINES.
C
C	A	(R) ARRAY A(INX,INY) CONTAINING VERTICAL HEIGHT DATA
C	INX,INY (I) DIMENSION OF A ARRAY
C	NX,NY	(I) INDICATING SIZE OF A ARRAY TO PLOT
C	PSI,PHI,THETA (R) ANGLES (IN DEGREES) OF INIT3DH
C	DV	(R) PERSPECTIVE FACTOR OF INIT3DH (9999=NO PERSPECTIVE)
C		    DV < 0 DO NOT INITIALIZE INIT3DH
C	XH,YH,ZH (R) LENGTH OF EACH AXIS (INCHES)
C	IFLAG	(I) PLOT FLAG
C		>1000	REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0	REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT)
C		= 0	CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0	SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C 	(ONE'S DIGIT) = 2 USE RAM TEK COLOR CONTROL ARRAY
C		      = 1 DO NOT USE RAM TEK COLOR ARRAY
C	(TEN'S DIGIT) = 0 NO HIDDEN LINE REMOVAL FOR AXIS
C		      = 1 HIDDEN LINE REMOVAL OF AXIS
C		      = 2 NO HIDDEN LINES AT ALL (SKETCH USED)
C		      = 3 NO HIDDEN LINES AT ALL (SKETCH NOT USED)
C 	(HUNDRED'S)   = 0 ASK WHICH SCREEN DEVICE TO USE
C		      <>0 SCREEN DEVICE CODE NUMBER
C 	IAXIS	(I) AXIS OPTION FLAG
C			= 0	DO NOT PLOT AXIS--FOLLOWING VARIABLES NOT REQUIRED
C			< 0	PLOT AXIS, USE INPUT Z AXIS SCALE
C			> 0	PLOT AXIS, USE COMPUTED Z AXIS SCALE
C 	(ONE'S DIGIT)	= 1	PLOT AXIS, Y AXIS SCALE--VARIABLES REQUIRED
C			= 2	PLOT AXIS, AUTO SCALE Y AXIS--VARIABLES REQUIRED
C 	XT,YT,ZT	(B) STRINGS FOR AXIS TITLES
C 	NXT,NYT,NZT	(I)  LENGTH OF AXIS TITLES
C			     IF ZERO THEN THAT AXIS NOT PLOTTED
C 	XASTART,YASTART	(R) AXIS START VALUES
C 	XAEND,YAEND	(R) AXIS END VALUES
C FOLLOWING ONLY REQUIRED IF IAXIS<>0
C	TX,TY,TZ	(R) NUMBER OF TICK MARKS (SEE AXIS3DH)
C	TSX,TSY,TSZ	(R) SIZE OF TITLE AND NUMBERS OF AXIS
C			  IF LESS THAN ZERO DO NOT AUTO-SCALE BY (x10^POWER)
C	PGX,PGY,PGZ	(R) ROTATION OF TEXT AND TICKS AROUND AXIS (DEG)
C	FSX,FSY,FSZ	(R) NUMBER FORMAT OF AXISES (SEE AXIS3DH)
C 	AMININ,AMAXIN 	(R) YAXIS SCALING FACTORS (ONLY NEEDED IF IAXIS < 0)
C	ICOL		(I) COLOR INDEX TABLE (REQUIRED IF 1'S DIGIT OF IFLAG=2
C				ICOL(1) AXIS COLOR
C				ICOL(2) SURFACE COLOR
C			    NOTE: WHEN BOTH AXIS AND SURFACE USE SKETCH
C			    BOTH ARE PLOTTED USING SURFACE COLOR
C
	DIMENSION A(INX,INY),PAS(2),ICOL(2)
	LOGICAL REPEAT/.FALSE./,HIDDEN
	INTEGER XT(1),YT(1),ZT(1)
C
	DIMENSION XX(4),YY(4),ZZ(4)
C
	PARAMETER (IWRKSZE=100000)
	COMMON /GO/ISIZE,WORK(IWRKSZE)
C
	IF (ISIZE.LT.IWRKSZE) ISIZE=IWRKSZE	! INCREASE WORKSPACE IF NEEDED
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	IF (IAXIS.LT.0) THEN
		AMAX=AMAXIN
		AMIN=AMININ
	ELSE
		AMAX=A(1,1)
		AMIN=A(1,1)
		DO 16 IY=1,NY			!DETERMINE MAX,MIN ARRAY VALUES
			DO 16 IX=1,NX
				AMAX=AMAX1(AMAX,A(IX,IY))
				AMIN=AMIN1(AMIN,A(IX,IY))
16		CONTINUE
	ENDIF
C
	XLEN=ABS(XH)
	YLEN=ABS(YH)
	ZLEN=ABS(ZH)
	IF (MOD(IABS(IAXIS),10).EQ.2) THEN
		PAS(1)=AMAX
		PAS(2)=AMIN
		CALL SCALE(PAS,ZLEN,2,1,1,AMIN,DAA)
		AMAX=ZLEN*DAA+AMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=MOD(JF,10000)/100
		CALL FRAME(3,ILU,1.5,0.65,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL PLOT MODE
	ENDIF
C
	HIDDEN=.TRUE.
	IOPT=MOD(JF,100)/10
	IF (IOPT.EQ.3) THEN
		IOPT=0
		HIDDEN=.FALSE.
	ENDIF
	JF=MOD(JF,10)
C
C	INITIALIZE 3DH PACKAGE
C
      	MNE=3			! MAX NUMBER OF EDGES
      	MNH=0			! NO HOLES
	ICORE=(25+5*MNE+4*MNH)*NX*NY
	IF (IAXIS.NE.0) ICORE=ICORE+10000	! ESTIMATED AXISES SPACE
	IF (ICORE.GT.ISIZE) THEN
		CALL CTERM(1)
		WRITE (*,3020)
3020		FORMAT(' *** WARNING: COMMON BLOCK /GO/ MAY BE TOO SMALL')
		CALL CTERM(-1)
	ENDIF
	SF=1.0			! SCALE FACTOR
	IF (IOPT.EQ.2) SF=-1.0	! NO HIDDEN LINE REMOVAL
	IF (DV.GT.0.) CALL INIT3DH(0.,0.,0.,PSI,PHI,THETA,SF,DV,MNE,MNH)
C
	IF (IAXIS.NE.0) THEN		!PLOT AXIS LABELS
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(1)),0.,0) !RAM TEK COLOR
		IF (NXT.GT.0) THEN 	!PLOT X AXIS
			CALL AXIS3DH(0.,0.,0.,0.,0.,PGX,XT,
     $				NXT,XLEN,XASTART,XAEND,
     $				TX,TSX,FSX,IOPT)
		ENDIF
		IF (NYT.GT.0) THEN 	!PLOT Y AXIS
			CALL AXIS3DH(0.,0.,0.,0.,90.,PGY,YT,
     $				NYT,YLEN,YASTART,YAEND,
     $				TY,TSY,FSY,IOPT)
		ENDIF
		IF (NZT.GT.0) THEN 	!PLOT Z AXIS
			CALL AXIS3DH(0.,0.,0.,90.,0.,PGZ,ZT,
     $				NZT,-ZLEN,AMIN,AMAX,
     $				TZ,TSZ,FSZ,IOPT)
		ENDIF
	ENDIF
	IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(2)),0.,0) !RAM TEK COLOR
C
C	BEGIN SURFACE PLOTTING
C
	XXX=NX
	DX=XLEN/IABS(NX)
	DY=YLEN/IABS(NY)
	DZ=1.0
	IF (AMAX-AMIN.NE.0.0) DZ=ZLEN/(AMAX-AMIN)
C
	DO 50 IY=1,IABS(NY)-1	! FOR EACH Y CELL
		Y1=(IY-1)*DY
		Y2=IY*DY
		DO 50 IX=1,IABS(NX)-1	! FOR EACH X CELL
			X1=(IX-1)*DX
			X2=IX*DX
C
			XX(1)=X1
			XX(2)=X2
			XX(3)=X2
			XX(4)=X1
			YY(1)=Y1
			YY(2)=Y1
			YY(3)=Y2
			YY(4)=Y1
			ZZ(1)=(A(IX  ,IY)   - AMIN)*DZ
			ZZ(2)=(A(IX+1,IY)   - AMIN)*DZ
			ZZ(3)=(A(IX+1,IY+1) - AMIN)*DZ
			ZZ(4)=ZZ(1)
			NC=4
			IERR=0
			IF (HIDDEN) THEN
				CALL SKETCH(XX,YY,ZZ,NC,IERR)
			ELSE
				IERR=3
				DO 20 JF=1,NC
				  CALL PLT3DH(XX(JF),YY(JF),ZZ(JF),IERR)
				  IERR=2
20				CONTINUE
			ENDIF
			IF (IERR.LT.0) THEN
				CALL CTERM(1)
		WRITE(*,3010) IERR
3010		FORMAT(' *** WARNING: SKETCH ERROR HAS OCCURED',I3)
				CALL CTERM(-1)
			ENDIF
C
			XX(1)=X1
			XX(2)=X1
			XX(3)=X2
			XX(4)=X1
			YY(1)=Y1
			YY(2)=Y2
			YY(3)=Y2
			YY(4)=Y1
			ZZ(1)=(A(IX  ,IY)   - AMIN)*DZ
			ZZ(2)=(A(IX  ,IY+1)   - AMIN)*DZ
			ZZ(3)=(A(IX+1,IY+1) - AMIN)*DZ
			ZZ(4)=ZZ(1)
			NC=4
			IERR=0
			IF (IY.EQ.IABS(NY)-1.AND.IX.EQ.IABS(NX)-1) IERR=1
			IF (HIDDEN) THEN
				CALL SKETCH(XX,YY,ZZ,NC,IERR)
			ELSE
				IERR=3
				DO 25 JF=1,NC
				  CALL PLT3DH(XX(JF),YY(JF),ZZ(JF),IERR)
				  IERR=2
25				CONTINUE
			ENDIF
			IF (IERR.LT.0) THEN
				CALL CTERM(1)
		WRITE (*,3200) IERR
3200		FORMAT(' *** WARNING: SKETCH ERROR HAS OCCURED',I3)
				CALL CTERM(-1)
			ENDIF
 50	CONTINUE
C
C	FINISH ROUTINE
C
	IF (.NOT.HIDDEN) CALL PLT3DH(XX(1),YY(1),ZZ(1),3)	!PEN UP
999	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL IN TEXT MODE
	ENDIF
	RETURN
	END
C
C ************************************************************************
C
	SUBROUTINE PLOTLGL(X,Y,W,NND,NPD,NN,NP,IFLAG,NSYM,
     1		XAXL,YAXL,NS,S,ILINE,NMX,NNX,MLX,TSX,NDX,
     2		SMX,NMY,NNY,MLY,TSY,NDY,SMY,XTITL,NXT,
     3		YTITL,NYT,TITLE,NT,XMIN,XMAX,YMIN,YMAX,ICOL)
C
C	CREATED BY D. LONG       MARCH, 1986	AT JPL
C
C	ROUTINE TO PLOT SEVERAL LINES OF DATA IN LOG/LINEAR COMBINATION
C	USING LINSEQ TO PROVIDE SOFTWARE LINE TYPING.  NOTE CAUTION IN
C	LINSEQ ROUTINE DESCRIPTION ABOUT CURVE FITTING WITH SPLINES.
C
C	X	REAL X VALUE ARRAY
C	Y	REAL Y DATA ARRAY DIMENSION Y(NPD,NND)
C	W	REAL WORKING ARRAY DIMENSIONED AT LEAST 3*NP+3
C	NND	INT  DIMENSION OF Y
C	NPD	INT  DIMENSION OF Y
C	NN	INT  # NUMBER OF LINES TO PLOT
C	NP	INT  # OF POINTS IN IN LINE
C	IFLAG	INT
C	  (MAG) >10000 REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0  REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT)
C		= 0  CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0  SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C	  (ONE'S DIGIT) = 1 X LINEAR, Y LOGRITHMIC (BASE 10) (SEE NOTE BELOW)
C	  (ONE'S DIGIT) = 2 X LOGRITHMIC,Y LINEAR
C	  (ONE'S DIGIT) = 3 X LOGRITHMIC,Y LOGRITHMIC
C	  (ONE'S DIGIT) = 4 X LINEAR, Y LINEAR
C	  (TEN'S DIGIT) = 0 NO AXIS OR GRID 
C	  (TEN'S DIGIT) = 1 PLOT BOX WITH AXIS TICKS
C	  (TEN'S DIGIT) = 2 PLOT SOLID GRID
C	  (TEN'S DIGIT) = 3 PLOT TICKED GRID W/O BOX
C	  (TEN'S DIGIT) = 4 PLOT TICKED GRID WITH BOX
C	  (TEN'S DIGIT) = 5 PLOT BOX WITH AXIS TICKS AND TICKED GRID
C	  (TEN'S DIGIT) = 6 PLOT W/O BOX OR GRID BUT USE AXISES
C	  (TEN'S DIGIT) = 7 PLOT SOLID LOGRITHMIC GRID
C	  (TEN'S DIGIT) = 8 PLOT DOTTED LOGRITHMIC GRID
C	  (TEN'S DIGIT) = 9 PLOT TICKED LOGRITHMIC GRID
C     (HUNDRED'S DIGIT) = 0 ASK WHICH SCREEN DEVICE TO USE
C			<>0 SCREEN DEVICE NUMBER CODE
C	NSYM	INT   NUMBER OF POINTS BETWEEN SYMBOLS
C		      < 0  SYMBOLS ONLY, NO LINE
C		      = 0  NO SYMBOLS, LINE ONLY
C		      > 0 NUMBER OF POINTS BETWEEN SYMBOLS
C	XAXL	REAL  X AXIS LENGTH IN INCHES
C		< 0  USE INPUT VALUE SCALING--AXIS PLOTTED
C		> 0  USE AUTO SCALING--AXIS PLOTTED
C	YAXL	REAL  Y AXIS LENGTH IN INCHES
C		< 0  USE INPUT VALUE SCALING--AXIS PLOTTED
C		> 0  USE AUTO SCALING--AXIS PLOTTED
C	NS	INT  NUMBER OF DATA SMOOTHING PASSES (NOMALLY ZERO)
C	S	REAL APPROXIMATE LINE TYPE PATTERN LENGTH (INCHES)
C	ILINE	INT  LINE TYPE ARRAY DIMENSIONED 5*NP
C			ILINE(1)= L1 FOR LINE 1
C			ILINE(2)= L2 FOR LINE 1
C			  ...        ...
C			ILINE(5)= L5 FOR LINE 1
C			ILINE(6)= L1 FOR LINE 2
C			ILINE(7)= L2 FOR LINE 3
C			  ...        ...	, ETC.
C
C NOTE: 	XMIN,XMAX,YMIN,YMAX ARE ONLY REQUIRED IF XAXL OR YAXL < 0
C
C	NMX...SMY    PARAMETERS FOR PLOT AXISES (SEE ROUTINE AXIS2)
C			ALL PARAMETERS MUST BE PRESENT AND WILL BE USED
C	XTITL	CHAR  X AXIS TITLE
C	NXT	INT   NUMBER OF CHARACTERS IN X AXIS TITLE STRING
C		< 0  TICKS ON INSIDE OF AXIS
C		= 0  NO AXIS
C		> 0  TICKS ON TITLE SIDE OF AXIS
C	YTITL	CHAR  Y AXIS TITLE
C	NYT	INT   NUMBER OF CHARACTERS IN Y AXIS TITLE STRING
C		< 0  TICKS ON INSIDE OF AXIS
C		= 0  NO AXIS
C		> 0  TICKS ON TITLE SIDE OF AXIS
C	TITLE	CHAR  PLOT TITLE
C	NT	INT   NUMBER OF CHARACTERS IN PLOT TITLE
C		= 0  NO TITLE
C		< 0  USE PEN COLOR ARRAY
C	XMIN	REAL  MINIMUM VALUE DISPLAYED ON X AXIS
C	XMAX	REAL  MAXIMUM VALUE DISPLAYED ON X AXIS
C	YMIN	REAL  MINIMUM VALUE OF Y ARRAY
C	YMAX	REAL  MAXIMUM VALUE OF Y ARRAY
C	ICOL	INT   ARRAY OF PEN COLORS PEN DIMENSIONED ICOL(5+NN)
C			ICOL(1) = GRID COLOR
C			ICOL(2) = AXIS LINE COLOR
C			ICOL(3) = AXIS NUMBERS COLOR
C			ICOL(4) = AXIS TITLE COLOR
C			ICOL(5) = AXIS EXPONENT COLOR
C			ICOL(6) = TITLE COLOR (COLOR UPON RETURN)
C			ICOL(7) = LINE COLOR 1
C			ICOL(8) = LINE COLOR 2
C			   .	     .	     3
C
C	NOTE: THE CONTENTS OF THE X AND Y ARRAYS ARE MODIFIED DURING CALL
C	RETURNED CONTENTS ARE THE PLOT SCALED VERSIONS OF INPUT VALUES
C
	INTEGER XTITL(1),YTITL(1),TITLE(1)
	REAL Y(NPD,NND),X(NPD),W(1),PY(2)
	LOGICAL REPEAT
	INTEGER PL,ICOL(1),IC(4),ILINE(1)
	COMMON /CPLOTLGL/XM,DX,YM,DY
	DATA REPEAT/.FALSE./
C
	PL=3					!PRINTER DATA FILE=FOR003.DAT
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND		!CLOSE PLOT PACKAGE
		REPEAT=.FALSE.
		RETURN
	ENDIF
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=-(JF/100)			!SCREEN DEVICE, CLEAR SCREEN
		CALL FRAME(PL,ILU,1.2,1.2,1.)	!INTIALIZE PLOT PACKAGE
	ELSE
		CALL CTERM(-1)			!RESTORE PLOT MODE
		JF=MOD(JF,100)
	ENDIF
	JF=MOD(JF,100)
	IF=JF
	IF (IF.GT.10) IF=MOD(IF,10)
	JF=JF-IF
	NX=ABS(ANINT(XAXL))				!GRID SIZES
	NY=ABS(ANINT(YAXL))
	IF (XAXL.GT.0.0) THEN				!INPUT SCALING
	    XM=X(1)
	    DX=X(1)
	    DO 5 I=2,NP
		XM=MAX(X(I),XM)
		DX=MIN(X(I),DX)
 5	    CONTINUE
	    PY(1)=XM
	    PY(2)=DX
	    IF (IF.EQ.2.OR.IF.EQ.3) THEN
		CALL SCALG(PY,ABS(XAXL),2,1,1,XM,DX)	!SMOOTH SCALE FACTORS
		GDX=1./DX
		NX=DX*ABS(XAXL)
		INX=-1
	    ELSE
		CALL SCALE(PY,ABS(XAXL),2,1,1,XM,DX)	!SMOOTH SCALE FACTORS
		GDX=1.
		INX=1
	    ENDIF
	ELSE
	    IF (IF.EQ.2.OR.IF.EQ.3) THEN
		XM=(ABS(XMIN)+1.E-35)
		IF (XM.NE.AINT(XM)) XM=AINT(XM)-1.
		DX=AINT(ALOG10(ABS(XMAX)+1.E-35))+1.
		DX=(DX-XM)/ABS(XAXL)
		GDX=1./DX
		NX=DX*ABS(XAXL)
		INX=-1
	    ELSE
		XM=XMIN
		DX=(XMAX-XMIN)/ABS(XAXL)
		GDX=1.
		INX=1
	    ENDIF
	ENDIF
	IF (YAXL.GT.0.0) THEN				!INPUT SCALING
	    YM=Y(1,1)
	    DY=Y(1,1)
	    DO 6 IN=1,NN
		DO 6 I=1,NP
		    YM=MAX(Y(I,IN),YM)
		    DY=MIN(Y(I,IN),DY)
6		CONTINUE
	    PY(1)=YM
	    PY(2)=DY
	    IF (IF.EQ.1.OR.IF.EQ.3) THEN
		CALL SCALG(PY,ABS(YAXL),2,1,1,YM,DY)	!SMOOTH SCALE FACTORS
		GDY=1./DY
		NY=DY*ABS(YAXL)
		INY=-1
	    ELSE
		CALL SCALE(PY,ABS(YAXL),2,1,1,YM,DY)	!SMOOTH SCALE FACTORS
		GDY=1.
		INY=1
	    ENDIF
	ELSE
	    IF (IF.EQ.1.OR.IF.EQ.3) THEN
		YM=ALOG10(ABS(YMIN)+1.E-35)
		IF (YM.NE.AINT(YM)) YM=AINT(YM)-1.
		DY=AINT(ALOG10(ABS(YMAX)+1.E-35))+1.
		DY=(DY-YM)/ABS(YAXL)
		GDY=1./DY
		NY=DY*ABS(YAXL)
		INY=-1
	    ELSE
		YM=YMIN
		DY=(YMAX-YMIN)/ABS(YAXL)
		GDY=1.
		INY=1
	    ENDIF
	ENDIF
	IF (JF.EQ.0) GOTO 77		!NO AXIS OR GRID
	IF (JF.EQ.60) GOTO 76		!NO GRID
	IF (JF.GT.60) THEN		!LOGRITHMIC GRID
		II=0
		IF (JF.EQ.70) II=1
		IF (JF.EQ.80) II=2
		CALL LGRID(0.,0.,GDX,GDY,INX*NX,INY*NY,II)
		GOTO 76
	ENDIF
	IF (JF.GE.20) THEN		!CARTESIAN GRID
		IF (JF.GE.30) NX=-NX
		IF (JF.GE.40) NY=-NY
		IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(1)),0.,0) !PEN COLOR
		CALL GRID(0.,0.,GDX,GDY,NX,NY)		!PRODUCE GRID
		IF (JF.EQ.50) JF=10
	ENDIF
76	NADD=0
	IF (NT.LT.0) THEN
		NADD=100000		!PEN COLOR
		IC(1)=ICOL(2)
		IC(2)=ICOL(3)
		IC(3)=ICOL(4)
		IC(4)=ICOL(5)
	ENDIF
	IX=IABS(NXT)
	XL=ABS(XAXL)*SIGN(1.,FLOAT(NXT))		!AXIS LENGTH WITH TICK INFOR
	IF (NXT.NE.0) THEN
		IF (IF.EQ.2.OR.IF.EQ.3) THEN
		CALL LGAXS(0.,0.,XTITL,-IX-NADD,XL,0.,XM,DX,IC) !X-AXIS
		IF (JF.EQ.10) CALL LGAXS(0.,ABS(YAXL),XTITL,IX+100+NADD,
     $			XL,0.,XM,DX,IC) !X-AXIS
		ELSE
		CALL AXIS2(0.,0.,XTITL,-IX-NADD-10000,XL,0.,XM,DX,
     $			NMX,NNX,-IABS(MLX),TSX,NDX,SMX,IC)  !X-AXIS
		IF (JF.EQ.10) CALL AXIS2(0.,ABS(YAXL),XTITL,IX+NADD+10100,
     $			XL,0.,XM,DX,NMX,NNX,-IABS(MLX),TSX,NDX,SMX,IC)  !X-AXIS
		ENDIF
	ENDIF
	IY=IABS(NYT)+1000
	YL=ABS(YAXL)*SIGN(1.,FLOAT(NYT))
	IF (NYT.NE.0) THEN
		IF (IF.EQ.1.OR.IF.EQ.3) THEN
		CALL LGAXS(0.,0.,YTITL,IY+NADD,YL,90.,YM,DY,IC) !Y-AXIS
		IF (JF.EQ.10) CALL LGAXS(ABS(XAXL),0.,YTITL,-IY-100-NADD,
     $			YL,90.,YM,DY,IC) !Y-AXIS
		ELSE
		CALL AXIS2(0.,0.,YTITL,IY+NADD+10000,YL,90.,YM,DY,
     $			NMY,NNY,-IABS(MLY),TSY,NDY,SMY,IC)  		!Y-AXIS
		IF (JF.EQ.10) CALL AXIS2(ABS(XAXL),0.,YTITL,-IY-NADD-10100,
     $			YL,90.,YM,DY,NMY,NNY,-IABS(MLY),TSY,NDY,SMY,IC)  !Y-AXIS
		ENDIF
	ENDIF
C
77	DO 78 I=1,NP
		IF (IF.EQ.2.OR.IF.EQ.3) X(I)=ALOG10(ABS(X(I))+1.E-35)
		X(I)=(X(I)-XM)/DX
78	CONTINUE
	DO 81 IN=1,NN
		DO 79 I=1,NP
			IF (IF.EQ.1.OR.IF.EQ.3) Y(I,IN)=ALOG10(ABS(Y(I,IN))
     $				+1.E-35)!LOG Y DATA
			Y(I,IN)=(Y(I,IN)-YM)/DY
79		CONTINUE
		IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(6+IN)),0.,0)!PEN COLOR
		I=(IN-1)*5+1
		IF (NSYM.GT.0) CALL LINSEQ(X,Y(1,IN),W,NP,NS,0.,S,ILINE(I),
     $			ILINE(I+1),ILINE(I+2),ILINE(I+3),ILINE(I+4))
		IF (NSYM.NE.0) THEN 
			DO 80 I=1,NP,IABS(NSYM)
				CALL SYMBOL(X(I),Y(I,IN),.2,IN,0.,1,-1)
80			CONTINUE
		ENDIF
81	CONTINUE
C
	CS=.21						!TITLE CHARACTER SIZE
	IF (NT.LT.0) CALL PLOT(FLOAT(ICOL(6)),0.,0)	!PEN COLOR
	IF (NT.NE.0.AND.JF.NE.0) THEN			!TITLE
		XP=0.0
		YP=0.0
		CALL SYMBOL(XP,YP,CS,TITLE,
     $			0.,MOD(IABS(NT),100),-3)	!TITLE LENGTH
		XP=ABS(XAXL/2.)-XP/2.
		IF (XP.LT.0.) XP=0.
		CALL SYMBOL(XP,ABS(YAXL)+.15,CS,TITLE,
     $			0.,MOD(IABS(NT),100),-1)	!PLOT TITLE
	ENDIF
	CALL PLOT(0.,0.,3)
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)				!ASK IF CLEAR SCREEN
		CALL PLOTND				!CLOSE PLOT PACKAGE
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)				!RETURN TO TEXT MODE
	ENDIF
	RETURN
	END
C
C **************************************************************************
C
	SUBROUTINE VAX5D(V,ND,N,A,B,IFLAG,IAX,IFORM,WORK,XH,YH,ZH,
     $		SPH,SPL,FAC,IW,IZ,ST,EN,T1,NT1,T2,NT2,T3,NT3,T4,NT4,
     $		VT,NVT,VM,VX,ICOL)
C
C THIS ROUTINE PLOTS A FOUR DIMENSIONAL FUNCTION AS A 2-D ARRAY OF 3-D SURFACES
C USING THE VAX3D SUBROUTINE.  COORDINATE SYSTEM USED IS SHOWN BELOW.  THE
C ROUTINE MAY PRODUCE SEVERAL PAGES OF OUTPUT WHICH CAN BE PASTED TOGETHER
C TO SEE THE ENTIRE ARRAY.  A FORMAT VALUE CAN BE USED TO INTERCHANGE THE
C DEPENDENT AXES OF THE PLOT TO PERMIT VIEWING THE FUNCTION IN DIFFERENT FORMS.
C A THREE DIMENSIONAL FUNCTION CAN BE PLOTTED BY USING A VALUE OF 1 FOR THE
C 4TH DIMENSION.
C
C COORDINATE SYSTEM USED:	IFORM=1 CORRESPONDS TO:
C
C       ^ W			  V(X,Y,Z,W)
C	|
C	|  V Y			EVERY POINT IN X,Y IS PLOTTED.
C	|  |/_X			ONLY THE IZth and IWth POINT (AS A NEW SURFACE)
C	|			Z,W IS PLOTTED. 
C	--------> Z
C
C	V	(R) ARRAY TO BE PLOTTED DIMENSIONED V(ND(1),ND(2),ND(3),ND(4))
C	ND	(I) ARRAY OF DIMENSIONS OF V
C	N	(I) ARRAY OF NUMBER OF POINTS IN EACH DIMENSION OF V TO PLOT
C	A,B	(R) ALPHA, BETA ANGLES FOR VAX3D
C	IFLAG	(I) OPTION FLAG
C			< 0 REPEAT PLOTTING, PLOT PACKAGE NOT CLOSED
C			= 0 PLOT PACKAGE CLOSED, NO PLOT PRODUCED
C			> 0 PLOT PACKAGE CLOSED, SINGLE PLOT PRODUCED
C	 (1000'S DIGIT) = 1 PLOT PACKAGE NOT INITIALIZED
C	 (ONE'S DIGIT)  = 2 COLOR INDEX ARRAY USED
C			= 1 COLOR TABLE IGNORED
C	 (TEN'S DIGIT)  = 0 PLOT SIDES
C			= 1 DO NOT PLOT SIDES
C	 (HUNDRED'S)	= 0 PROMPT FOR SCREEN GRAPHICS TYPE
C			<>0 SCREEN DEVICE NUMBER (SEE FRAME)
C	IAX	(I) AXIS OPTION FLAG
C			< 0 PLOT AXIS USING INPUT SCALE FACTORS VM,VX
C			= 0 AXIS NOT PLOTTED
C			> 0 PLOT AXIS USING MAX, MIN OF V ARRAY
C	 (ONE'S DIGIT)	= 1 USE ACTUAL MAX, MIN VALUES
C			= 2 USED SMOOTHED MAX, MIN VALUES
C	 (TEN'S DIGIT)  = 0 AXES PLOTTED FOR ALL SURFACES
C			= 1 AXES PLOTTED ONLY FOR FIRST SURFACE
C	IFORM	(I) PLOT FORMAT CODE
C			    AXIS USED FOR:	
C		     CODE     X  Y  Z  W
C			 1 :  1  2  3  4
C			 2 :  1  2  4  3
C			 3 :  1  4  2  3
C			 4 :  1  4  3  2
C			 5 :  1  3  2  4
C			 6 :  1  3  4  2
C			 7 :  2  1  3  4
C			 8 :  2  1  4  3
C			 9 :  2  4  3  1
C			10 :  2  4  1  3
C			11 :  2  3  4  1
C			12 :  2  3  1  4
C			13 :  3  1  2  4
C			14 :  3  1  4  2
C			15 :  3  2  1  4
C			16 :  3  2  4  1
C			17 :  3  4  1  2
C			18 :  3  4  2  1
C			19 :  4  1  2  3
C			20 :  4  1  3  2
C			21 :  4  2  1  3
C			22 :  4  2  3  1
C			23 :  4  3  1  2
C			24 :  4  3  2  1
C	WORK	(R) WORKING ARRAY DIMENSIONED N(X)*N(Y)
C	XH,YH,ZH (R) PLOTTING DIMENSIONS FOR SURFACE IN VAX3D
C	SPH,SPL	(R) SIZE OF PLOT PAGE (HEIGHT,LENGTH)
C	FAC	(R) SHRINK FACTOR FOR SURFACES (2=FACTOR(1/2))
C	IW,IZ	(R) PLOT EVERY IWth and IZth SURFACE
C	ST	(R) ARRAY OF STARTING VALUES FOR AXES FOR EACH DIMENSION OF V
C	EN	(R) ARRAY OF ENDING VALUES FOR AXES FOR EACH DIMENSION OF V
C	T1,T2,T3,T4,VT (B) TITLES FOR EACH AXIS OF EACH DIMENSION OF V
C	NT1,NT2,NT3,NT4,NVT (B) NUMBER OF CHARACTERS IN EACH TITLE
C OPTIONAL PARAMETERS REQUIRED ONLY IF IAX < 0 OR IFLAG ONE'S DIGIT=2
C	VM,VX	(R) MINIMUM AND MAXIMUM OF PLOT
C	ICOL	(I) ARRAY OF COLOR TABLE INDEXES
C
	DIMENSION ND(4)
	DIMENSION V(ND(1),ND(2),ND(3),ND(4)),N(4),
     $		WORK(1),ST(4),EN(4),ICOL(1)
	INTEGER T1(1),T2(1),T3(1),T4(1),VT(1)
C
	INTEGER I(4),XT,ZT,WT,TT
	INTEGER X(24),Y(24),Z(24),W(24)
	DATA X/1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4/
	DATA Y/2,2,3,3,4,4,1,1,4,4,3,3,1,1,2,2,4,4,1,1,2,2,3,3/
	DATA Z/3,4,2,4,2,4,3,4,3,1,4,1,2,4,1,4,1,2,2,3,1,3,1,2/
	DATA W/4,3,3,2,4,2,4,3,1,3,1,4,4,2,4,1,2,1,3,2,3,1,2,1/
C
	DIMENSION AS(2)
	LOGICAL REPEAT
	DATA REPEAT/.FALSE./,DTR/0.0174532925/
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	IF (IAX.LT.0) THEN	! USE INPUT MAX,MIN
		AMAX=VX
		AMIN=VM
	ELSE			! DETERMINE MAX,MIN ARRAY VALUES
		AMAX=V(1,1,1,1)
		AMIN=AMAX
		DO 16 I1=1,N(1)
			DO 16 I2=1,N(2)
				DO 16 I3=1,N(3)
					DO 16 I4=1,N(4)
				AMAX=AMAX1(AMAX,V(I1,I2,I3,I4))
				AMIN=AMIN1(AMIN,V(I1,I2,I3,I4))
16		CONTINUE
	ENDIF
	IF (MOD(IABS(IAXIS),10).EQ.2) THEN	! SMOOTH SCALE FACTORS
		AS(1)=AMAX
		AS(2)=AMIN
		CALL SCALE(AS,YH,2,1,1,AMIN,DAA)
		AMAX=YLEN*DAA+AMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=JF/100
		CALL FRAME(3,ILU,.8,.4,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL IN PLOT MODE
	ENDIF
C
	IF (IFORM.LT.1.OR.IFORM.GT.24) GOTO 1000
	JFF=-(10000+MOD(JF,100))
	IA=-1
	SP=0.3						!SPACING BETWEEN SURFACES
	IF (IAX.NE.0.AND.IAX/10.EQ.0) SP=1.2
C
C	REARRANGE AXISES.  THIS TECHNIQUE IS VERY MACHINE DEPENDENT.
C
C	%LOC() RETURNS THE ADDRESS OF THE VARIABLE
C
	I(1)=%LOC(T1(1))
	I(2)=%LOC(T2(1))
	I(3)=%LOC(T3(1))
	I(4)=%LOC(T4(1))
	XT=I(X(IFORM))
	ZT=I(Y(IFORM))
	WT=I(W(IFORM))
	TT=I(Z(IFORM))
	I(1)=NT1
	I(2)=NT2
	I(3)=NT3
	I(4)=NT4
	NNX=I(X(IFORM))
	NNZ=I(Y(IFORM))
	NNW=I(W(IFORM))
	NNT=I(Z(IFORM))
	NNV=NVT
	IF (IAX.EQ.0) THEN
		NNX=0.
		NNZ=0.
		NNV=0.
	ENDIF
C
	XS=ST(X(IFORM))
	XE=EN(X(IFORM))
	YS=ST(Y(IFORM))
	YE=EN(Y(IFORM))
C
C	DETERMINE THE NUMBER OF PAGES
C
	NX=N(X(IFORM))
	NY=N(Y(IFORM))
	NZ=(N(Z(IFORM))-1)/IZ+1
	NW=(N(W(IFORM))-1)/IW+1
	NZ1=NZ
	IF (NZ1.GT.1) NZ1=NZ-1
	NW1=NW
	IF (NW1.GT.1) NW1=NW-1
C
	AL=XH*COS(A*DTR)+COS(B*DTR)*ZH+SP	! SURFACE LENGTH
	AH=XH*SIN(A*DTR)+ZH*SIN(B*DTR)+YH+SP	! SURFACE HEIGHT
	XX=XH*COS(A*DTR)+0.5
	NPZ=IFIX((AL*NZ)/SPL/FAC)+1
	NPW=IFIX((AH*NW)/SPH/FAC)+1
C
C	SET SCALE FACTOR DOWN
C
	CALL FACTOR(1./FAC)
C
C	MAKE A PASS FOR EACH PAGE
C
	DO 500 IPW=1,NPW
	  IF (IPW.GT.1) THEN
		CALL PLOT(0.,0.,3)
		CALL NEWPAGE
		CALL CTERM(2)
		CALL RTERM(2)
		CALL CTERM(-1)
	  ENDIF
	  XWP=SPH*(IPW-1)*FAC
	DO 500 IPZ=1,NPZ
	  IF (IPZ.GT.1) THEN
		CALL PLOT(0.,0.,3)
		CALL NEWPAGE
		CALL CTERM(2)
		CALL RTERM(2)
		CALL CTERM(-1)
	  ENDIF
	  YZP=SPL*(IPZ-1)*FAC
	IF (IPZ.EQ.1.AND.IPW.EQ.1) THEN
	  CALL SYMBOL(-.49,0.0,.2,%VAL(WT),90.,NNW,-1)
	  CALL SYMBOL(0.2,-.3,.2,%VAL(TT),0.,NNT,-1)
	  CALL SYMBOL(999.,999.,.2,'  Page (',0.,8,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(IPZ),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(IPX),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,') of (',0.,6,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(NPZ),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(NPX),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,')',0.,1,-1)
	ENDIF
C
C	FOR EACH W AND Z POINT, CREATE A 3-D PLOT
C
	DO 300 KW=1,NW
	  JW=1+(KW-1)*IW
	  I(W(IFORM))=JW
	DO 300 KZ=1,NZ
	  JZ=1+(KZ-1)*IZ
	  I(Z(IFORM))=JZ
C
	YZ=AL*(KZ-1)-YZP
	XW=AH*(KW-1)-XWP
C
C	CHECK TO SEE IF WE DON'T NEED TO PLOT THIS SURFACE
C
	IF (YZ.GT.FAC*SPL.OR.YZ.LT.-AL) GOTO 300
	IF (XW.GT.FAC*SPH.OR.XW.LT.-AH) GOTO 300
C
C	MOVE PLOT ORIGIN
C
	CALL PLOT(YZ,XW,-3)
C
C	LABEL PLOT
C
	CALL CTERM(-1)
	CALL SYMBOL(XX,0.,.2,'(',0.,1,-1)
	CALL NUMBER(999.,999.,.2,FLOAT(JZ),0.,-1,-1)
	CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	CALL NUMBER(999.,999.,.2,FLOAT(JW),0.,-1,-1)
	IF (IAX.NE.0) THEN
		WVAL=(JW-1)*(EN(W(IFORM))-ST(W(IFORM)))/NW1+ST(W(IFORM))
		ZVAL=(JZ-1)*(EN(Z(IFORM))-ST(Z(IFORM)))/NZ1+ST(Z(IFORM))
		CALL SYMBOL(999.,999.,.2,') (',0.,3,-1)
		CALL NUMBER(999.,999.,.2,ZVAL,0.,2,-1)
		CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
		CALL NUMBER(999.,999.,.2,WVAL,0.,2,-1)
	ENDIF
	CALL SYMBOL(999.,999.,.2,')',0.,1,-1)
C
C	FILL WORKING ARRAY
C
	DO 150 IY=1,NY
		I(Y(IFORM))=IY
		DO 150 IX=1,NX
			I(X(IFORM))=IX
			WORK(IY+NY*(IX-1))=V(I(1),I(2),I(3),I(4))
150	CONTINUE
C
C	CALL VAX3D
C
	IF (MOD(JFF,10).EQ.2) THEN
	   CALL VAX3D(WORK,NX,NY,NX,NY,A,B,XH,YH,ZH,JFF,IA,%VAL(XT),NNX,
     $		XS,XE,VT,NNV,%VAL(ZT),NNZ,YS,YE,AMIN,AMAX,ICOL)
	ELSE
	   CALL VAX3D(WORK,NX,NY,NX,NY,A,B,XH,YH,ZH,JFF,IA,%VAL(XT),NNX,
     $		XS,XE,VT,NNV,%VAL(ZT),NNZ,YS,YE,AMIN,AMAX)
	ENDIF
C
C	TURN OFF AXES AFTER FIRST PLOT IF DESIRED
C
	IF (IAX/10.NE.0) THEN
		NNX=0.
		NNZ=0.
		NNV=0.
	ENDIF
C
C	RESTORE PLOT ORIGIN
C
	CALL PLOT(-YZ,-XW,-3)
300	CONTINUE
500	CONTINUE
C
C	RESTORE SCALE FACTOR
C
	CALL FACTOR(FAC)
C
C	FINISH UP PLOT
C
1000	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL IN TEXT MODE
	ENDIF
	RETURN
	END
C
C **************************************************************************
C
	SUBROUTINE MVAX5D(V,ND,N,A,B,IFLAG,IAX,IFORM,WORK,XH,YH,ZH,
     $		SPH,SPL,FAC,IW,IZ,ST,EN,T1,NT1,T2,NT2,T3,NT3,T4,NT4,
     $		VT,NVT,VM,VX,ICOL)
C
C THIS ROUTINE PLOTS A FOUR DIMENSIONAL FUNCTION AS A 2-D ARRAY OF 3-D SURFACES
C USING THE MVAX3D SUBROUTINE.  COORDINATE SYSTEM USED IS SHOWN BELOW.  THE
C ROUTINE MAY PRODUCE SEVERAL PAGES OF OUTPUT WHICH CAN BE PASTED TOGETHER
C TO SEE THE ENTIRE ARRAY.  A FORMAT VALUE CAN BE USED TO INTERCHANGE THE
C DEPENDENT AXES OF THE PLOT TO PERMIT VIEWING THE FUNCTION IN DIFFERENT FORMS.
C A THREE DIMENSIONAL FUNCTION CAN BE PLOTTED BY USING A VALUE OF 1 FOR THE
C 4TH DIMENSION.
C
C COORDINATE SYSTEM USED:	IFORM=1 CORRESPONDS TO:
C
C       ^ W			  V(X,Y,Z,W)
C	|
C	|   V 			EVERY POINT IN X,Y IS PLOTTED.
C	|   |_X			ONLY THE IZth and IWth POINT (AS A NEW SURFACE)
C	| Y/			Z,W IS PLOTTED. 
C	--------> Z
C
C	V	(R) ARRAY TO BE PLOTTED DIMENSIONED V(ND(1),ND(2),ND(3),ND(4))
C	ND	(I) ARRAY OF DIMENSIONS OF V
C	N	(I) ARRAY OF NUMBER OF POINTS IN EACH DIMENSION OF V TO PLOT
C	A,B	(R) ALPHA, BETA ANGLES FOR MVAX3D
C	IFLAG	(I) OPTION FLAG
C			< 0 REPEAT PLOTTING, PLOT PACKAGE NOT CLOSED
C			= 0 PLOT PACKAGE CLOSED, NO PLOT PRODUCED
C			> 0 PLOT PACKAGE CLOSED, SINGLE PLOT PRODUCED
C	 (1000'S DIGIT) = 1 PLOT PACKAGE NOT INITIALIZED
C	 (ONE'S DIGIT)  = 2 COLOR INDEX ARRAY USED
C			= 1 COLOR TABLE IGNORED
C	 (TEN'S DIGIT)  = 0 PLOT MESH WITH SIDES, LOWER SURFACE NOT SHOWN
C			= 1 MESH W/O SIDES, LOWER SURFACE DOTTED
C			= 2 MESH W/O SIDES, LOWER SURFACE SOLID
C			= 3 HISTOGRAM COLUMS, LOWER SURFACE DOTTED
C			= 4 HISTOGRAM COLUMS, LOWER SURFACE SOLID
C	 (HUNDRED'S)	= 0 PROMPT FOR SCREEN GRAPHICS TYPE
C			<>0 SCREEN DEVICE NUMBER (SEE FRAME)
C	IAX	(I) AXIS OPTION FLAG
C			< 0 PLOT AXIS USING INPUT SCALE FACTORS VM,VX
C			= 0 AXIS NOT PLOTTED
C			> 0 PLOT AXIS USING MAX, MIN OF V ARRAY
C	 (ONE'S DIGIT)	= 1 USE ACTUAL MAX, MIN VALUES
C			= 2 USED SMOOTHED MAX, MIN VALUES
C	 (TEN'S DIGIT)  = 0 AXES PLOTTED FOR ALL SURFACES
C			= 1 AXES PLOTTED ONLY FOR FIRST SURFACE
C	IFORM	(I) PLOT FORMAT CODE
C			    AXIS USED FOR:	
C		     CODE     X  Y  Z  W
C			 1 :  1  2  3  4
C			 2 :  1  2  4  3
C			 3 :  1  4  2  3
C			 4 :  1  4  3  2
C			 5 :  1  3  2  4
C			 6 :  1  3  4  2
C			 7 :  2  1  3  4
C			 8 :  2  1  4  3
C			 9 :  2  4  3  1
C			10 :  2  4  1  3
C			11 :  2  3  4  1
C			12 :  2  3  1  4
C			13 :  3  1  2  4
C			14 :  3  1  4  2
C			15 :  3  2  1  4
C			16 :  3  2  4  1
C			17 :  3  4  1  2
C			18 :  3  4  2  1
C			19 :  4  1  2  3
C			20 :  4  1  3  2
C			21 :  4  2  1  3
C			22 :  4  2  3  1
C			23 :  4  3  1  2
C			24 :  4  3  2  1
C	WORK	(R) WORKING ARRAY DIMENSIONED N(X)*N(Y)
C	XH,YH,ZH (R) PLOTTING DIMENSIONS FOR SURFACE IN MVAX3D
C	SPH,SPL	(R) SIZE OF PLOT PAGE (HEIGHT,LENGTH)
C	FAC	(R) SHRINK FACTOR FOR SURFACES (2=FACTOR(1/2))
C	IW,IZ	(R) PLOT EVERY IWth and IZth SURFACE
C	ST	(R) ARRAY OF STARTING VALUES FOR AXES FOR EACH DIMENSION OF V
C	EN	(R) ARRAY OF ENDING VALUES FOR AXES FOR EACH DIMENSION OF V
C	T1,T2,T3,T4,VT (B) TITLES FOR EACH AXIS OF EACH DIMENSION OF V
C	NT1,NT2,NT3,NT4,NVT (B) NUMBER OF CHARACTERS IN EACH TITLE
C OPTIONAL PARAMETERS REQUIRED ONLY IF IAX < 0 OR IFLAG ONE'S DIGIT=2
C	VM,VX	(R) MINIMUM AND MAXIMUM OF PLOT
C	ICOL	(I) ARRAY OF COLOR TABLE INDEXES
C
	DIMENSION ND(4)
	DIMENSION V(ND(1),ND(2),ND(3),ND(4)),N(4),
     $		WORK(1),ST(4),EN(4),ICOL(1)
	INTEGER T1(1),T2(1),T3(1),T4(1),VT(1)
C
	INTEGER I(4),XT,ZT,WT,TT
	INTEGER X(24),Y(24),Z(24),W(24)
	DATA X/1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4/
	DATA Y/2,2,3,3,4,4,1,1,4,4,3,3,1,1,2,2,4,4,1,1,2,2,3,3/
	DATA Z/3,4,2,4,2,4,3,4,3,1,4,1,2,4,1,4,1,2,2,3,1,3,1,2/
	DATA W/4,3,3,2,4,2,4,3,1,3,1,4,4,2,4,1,2,1,3,2,3,1,2,1/
C
	DIMENSION AS(2)
	LOGICAL REPEAT
	DATA REPEAT/.FALSE./,DTR/0.0174532925/
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND	!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	IF (IAX.LT.0) THEN	! USE INPUT MAX,MIN
		AMAX=VX
		AMIN=VM
	ELSE			! DETERMINE MAX,MIN ARRAY VALUES
		AMAX=V(1,1,1,1)
		AMIN=AMAX
		DO 16 I1=1,N(1)
			DO 16 I2=1,N(2)
				DO 16 I3=1,N(3)
					DO 16 I4=1,N(4)
				AMAX=AMAX1(AMAX,V(I1,I2,I3,I4))
				AMIN=AMIN1(AMIN,V(I1,I2,I3,I4))
16		CONTINUE
	ENDIF
	IF (MOD(IABS(IAXIS),10).EQ.2) THEN	! SMOOTH SCALE FACTORS
		AS(1)=AMAX
		AS(2)=AMIN
		CALL SCALE(AS,YH,2,1,1,AMIN,DAA)
		AMAX=YLEN*DAA+AMIN
	ENDIF
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=JF/100
		CALL FRAME(3,ILU,.8,.4,1.)		!INTIALIZE
	ELSE
		CALL CTERM(-1)				!PUT TERMINAL IN PLOT MODE
	ENDIF
C
	IF (IFORM.LT.1.OR.IFORM.GT.24) GOTO 1000
	JFF=-(10000+MOD(JF,100))
	IA=-1
	SP=0.3						!SPACING BETWEEN SURFACES
	IF (IAX.NE.0.AND.IAX/10.EQ.0) SP=1.2
C
C	REARRANGE AXISES.  THIS TECHNIQUE IS VERY MACHINE DEPENDENT.
C
C	%LOC() RETURNS THE ADDRESS OF THE VARIABLE
C
	I(1)=%LOC(T1(1))
	I(2)=%LOC(T2(1))
	I(3)=%LOC(T3(1))
	I(4)=%LOC(T4(1))
	XT=I(X(IFORM))
	ZT=I(Y(IFORM))
	WT=I(W(IFORM))
	TT=I(Z(IFORM))
	I(1)=NT1
	I(2)=NT2
	I(3)=NT3
	I(4)=NT4
	NNX=I(X(IFORM))
	NNZ=I(Y(IFORM))
	NNW=I(W(IFORM))
	NNT=I(Z(IFORM))
	NNV=NVT
	IF (IAX.EQ.0) THEN
		NNX=0.
		NNZ=0.
		NNV=0.
	ENDIF
C
	XS=ST(X(IFORM))
	XE=EN(X(IFORM))
	YS=ST(Y(IFORM))
	YE=EN(Y(IFORM))
C
C	DETERMINE THE NUMBER OF PAGES
C
	NX=N(X(IFORM))
	NY=N(Y(IFORM))
	NZ=(N(Z(IFORM))-1)/IZ+1
	NW=(N(W(IFORM))-1)/IW+1
	NZ1=NZ
	IF (NZ1.GT.1) NZ1=NZ-1
	NW1=NW
	IF (NW1.GT.1) NW1=NW-1
C
	AL=XH*COS(A*DTR)+COS(B*DTR)*ZH+SP	! SURFACE LENGTH
	AH=XH*SIN(A*DTR)+ZH*SIN(B*DTR)+YH+SP	! SURFACE HEIGHT
	XX=XH*COS(A*DTR)+0.5
	NPZ=IFIX((AL*NZ)/SPL/FAC)+1
	NPW=IFIX((AH*NW)/SPH/FAC)+1
C
C	SET SCALE FACTOR DOWN
C
	CALL FACTOR(1./FAC)
C
C	MAKE A PASS FOR EACH PAGE
C
	DO 500 IPW=1,NPW
	  IF (IPW.GT.1) THEN
		CALL PLOT(0.,0.,3)
		CALL NEWPAGE
		CALL CTERM(2)
		CALL RTERM(2)
		CALL CTERM(-1)
	  ENDIF
	  XWP=SPH*(IPW-1)*FAC
	DO 500 IPZ=1,NPZ
	  IF (IPZ.GT.1) THEN
		CALL PLOT(0.,0.,3)
		CALL NEWPAGE
		CALL CTERM(2)
		CALL RTERM(2)
		CALL CTERM(-1)
	  ENDIF
	  YZP=SPL*(IPZ-1)*FAC
	IF (IPZ.EQ.1.AND.IPW.EQ.1) THEN
	  CALL SYMBOL(-.49,0.0,.2,%VAL(WT),90.,NNW,-1)
	  CALL SYMBOL(0.2,-.3,.2,%VAL(TT),0.,NNT,-1)
	  CALL SYMBOL(999.,999.,.2,'  Page (',0.,8,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(IPZ),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(IPX),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,') of (',0.,6,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(NPZ),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	  CALL NUMBER(999.,999.,.2,FLOAT(NPX),0.,-1,-1)
	  CALL SYMBOL(999.,999.,.2,')',0.,1,-1)
	ENDIF
C
C	FOR EACH W AND Z POINT, CREATE A 3-D PLOT
C
	DO 300 KW=1,NW
	  JW=1+(KW-1)*IW
	  I(W(IFORM))=JW
	DO 300 KZ=1,NZ
	  JZ=1+(KZ-1)*IZ
	  I(Z(IFORM))=JZ
C
	YZ=AL*(KZ-1)-YZP
	XW=AH*(KW-1)-XWP
C
C	CHECK TO SEE IF WE DON'T NEED TO PLOT THIS SURFACE
C
	IF (YZ.GT.FAC*SPL.OR.YZ.LT.-AL) GOTO 300
	IF (XW.GT.FAC*SPH.OR.XW.LT.-AH) GOTO 300
C
C	MOVE PLOT ORIGIN
C
	CALL PLOT(YZ,XW,-3)
C
C	LABEL PLOT
C
	CALL CTERM(-1)
	CALL SYMBOL(XX,0.,.2,'(',0.,1,-1)
	CALL NUMBER(999.,999.,.2,FLOAT(JZ),0.,-1,-1)
	CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
	CALL NUMBER(999.,999.,.2,FLOAT(JW),0.,-1,-1)
	IF (IAX.NE.0) THEN
		WVAL=(JW-1)*(EN(W(IFORM))-ST(W(IFORM)))/NW1+ST(W(IFORM))
		ZVAL=(JZ-1)*(EN(Z(IFORM))-ST(Z(IFORM)))/NZ1+ST(Z(IFORM))
		CALL SYMBOL(999.,999.,.2,') (',0.,3,-1)
		CALL NUMBER(999.,999.,.2,ZVAL,0.,2,-1)
		CALL SYMBOL(999.,999.,.2,',',0.,1,-1)
		CALL NUMBER(999.,999.,.2,WVAL,0.,2,-1)
	ENDIF
	CALL SYMBOL(999.,999.,.2,')',0.,1,-1)
C
C	FILL WORKING ARRAY
C
	DO 150 IY=1,NY
		I(Y(IFORM))=IY
		DO 150 IX=1,NX
			I(X(IFORM))=IX
			WORK(IY+NY*(IX-1))=V(I(1),I(2),I(3),I(4))
150	CONTINUE
C
C	CALL MVAX3D
C
	IF (MOD(JFF,10).EQ.2) THEN
	   CALL MVAX3D(WORK,NX,NY,NX,NY,A,B,XH,YH,ZH,JFF,IA,%VAL(XT),
     $		NNX,XS,XE,VT,NNV,%VAL(ZT),NNZ,YS,YE,AMIN,AMAX,ICOL)
	ELSE
	   CALL MVAX3D(WORK,NX,NY,NX,NY,A,B,XH,YH,ZH,JFF,IA,%VAL(XT),
     $		NNX,XS,XE,VT,NNV,%VAL(ZT),NNZ,YS,YE,AMIN,AMAX)
	ENDIF
C
C	TURN OFF AXES AFTER FIRST PLOT IF DESIRED
C
	IF (IAX/10.NE.0) THEN
		NNX=0.
		NNZ=0.
		NNV=0.
	ENDIF
C
C	RESTORE PLOT ORIGIN
C
	CALL PLOT(-YZ,-XW,-3)
300	CONTINUE
500	CONTINUE
C
C	RESTORE SCALE FACTOR
C
	CALL FACTOR(FAC)
C
C	FINISH UP PLOT
C
1000	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL IN TEXT MODE
	ENDIF
	RETURN
	END
C
C *************************************************************************
C
	SUBROUTINE CVAX3D(A,INX,INZ,NX,NZ,ALPHA,BETA,XH,YH,ZH,IFLAG,
     $		IAXIS,CLOW,NLEV,CVAL,NCH,NMIN,ILINE,IPLS,
     $		XT,NXT,XASTART,XAEND,YT,NYT,ZT,NZT,ZASTART,ZAEND,
     $		AMININ,AMAXIN,ICOL)
C
C	CREATED BY D. LONG  AUG 1987	AT JPL
C
C	SIMPLIFIES CALLING CVAX3DX ROUTINE
C
	DIMENSION A(INX,INZ),ICOL(1),ILINE(1),CVAL(1),IPLS(1)
	INTEGER XT(1),YT(1),ZT(1)
	IAX=MOD(IAXIS,100)
	IF (IAXIS.LT.0) THEN
		AMN=AMININ
		AMX=AMAXIN
	ENDIF
	IF (MOD(IABS(IFLAG),10).EQ.2) THEN
		CALL CVAX3DX(A,INX,INZ,NX,NZ,ALPHA,BETA,XH,YH,ZH,IFLAG,
     $		        IAX,CLOW,NLEV,CVAL,NCH,NMIN,ILINE,IPLS,
     $			XT,NXT,XASTART,XAEND,TX,TSX,FDX,
     $			YT,NYT,TY,TSY,FDY,
     $			ZT,NZT,ZASTART,ZAEND,TZ,TSZ,FDZ,
     $			AMN,AMX,ICOL)
	ELSE
		CALL CVAX3DX(A,INX,INZ,NX,NZ,ALPHA,BETA,XH,YH,ZH,IFLAG,
     $		        IAX,CLOW,NLEV,CVAL,NCH,NMIN,ILINE,IPLS,
     $			XT,NXT,XASTART,XAEND,TX,TSX,FDX,
     $			YT,NYT,TY,TSY,FDY,
     $			ZT,NZT,ZASTART,ZAEND,TZ,TSZ,FDZ,
     $			AMN,AMX,IC)
	ENDIF
	RETURN
	END
C
C *********************************************************************
C
	SUBROUTINE CVAX3DX(A,INX,INZ,NX,NZ,ALPHA,BETA,XH,YH,ZH,IFLAG,
     $		IAXIS,CLOW,NLEV,CVAL,NCH,NMIN,ILINE,IPLS,
     $		XT,NXT,XASTART,XAEND,TX,TSX,FDX,
     $		YT,NYT,TY,TSY,FDY,
     $		ZT,NZT,ZASTART,ZAEND,TZ,TSZ,FDZ,
     $		AMININ,AMAXIN,ICOL)
C
C	CREATED BY D. LONG     AUG 1987	AT JPL
C
C	ROUTINE TO PLOT DATA IN 3-D MESH OR HISTOGRAM COLUMN FORM
C	WITH HIDDEN LINE REMOVAL USING PLT3D OR HLT3D
C	WITH A CONTOUR PLOT UNDERNEATH USING GCONTR
C
C	COORDINATE SYSTEM IS:		Y
C			VARIABLE NAMES:	|	
C				      Z/ \X
C
C	A	REAL ARRAY A(INX,INZ) CONTAINING VERTICAL HEIGHT DATA
C	INX,INZ INTEGERS DIMENSION OF A ARRAY
C	NX,NZ	INTEGERS INDICATING SIZE OF A ARRAY TO PLOT
C	ALPHA	REAL ANGLE (IN DEGREES) OF X AXIS (NX) (ALTITUDE)
C	BETA	REAL ANGLE (IN DEGREES) OF Z AXIS (NZ) (AZIMUTH)
C	XH,YH,ZH REAL LENGTH OF EACH AXIS
C	IFLAG	INTEGER
C		>10000	REPEAT PLOTTING BUT DO NOT INTIALIZE PLOTTER
C		< 0	REPEAT PLOTTING (DON'T CLOSE PLOTTER OUTPUT
C		= 0	CLOSE PLOTTER OUTPUT--DO NOT PRODUCE PLOT
C		> 0	SINGLE PLOT OUTPUT ONLY (PLOTTER OUTPUT CLOSED)
C 	(ONE'S DIGIT) = 2 USE PEN COLOR CONTROL ARRAY
C		      = 1 DO NOT USE PEN COLOR ARRAY
C	(TEN'S)       = 0 MESH (PLT3D) USED
C		      = 1 HISTOGRAM (MLT3D) USED
C 	(HUNDREDS'S DIGIT) = 0 ASK WHICH SCREEN DEVICE TO USE
C		           <>0 SCREEN DEVICE CODE NUMBER
C 	IAXIS	INTEGER AXIS OPTION FLAG
C			= 0	DO NOT PLOT AXIS--FOLLOWING VARIABLES NOT REQUIRED
C			< 0	PLOT AXIS, USE INPUT Y AXIS SCALE--FOLLOWING VARIABLES REQUIRED
C			> 0	PLOT AXIS, USE COMPUTED Y AXIS SCALE--FOLLOWING VARIABLES REQUIRED
C 	(ONE'S DIGIT)	= 1	PLOT AXIS, Y AXIS SCALE--VARIABLES REQUIRED
C			= 2	PLOT AXIS, AUTO SCALE Y AXIS--VARIABLES REQUIRED
C	(TEN'S DIGIT)	= 0	CONTOUR AXES, MESH AXES, BACK PLANE PLOTTED
C			= 1	CONTOUR AXES, MESH AXES PLOTTED W/O BACK PLANE
C			= 2	CONTOUR AXES PLOTTED W/O MESH AXES 
C			= 3	MESH AXES PLOTTED W/O BACK PLANE OR CONTOUR AXES
C			= 4	MESH AXES AND BACK PLANE W/O CONTOUR AXES
C	(HUNDRED'S)	= 0	DEFAULT AXIS PARAMETERS
C			= 1	SPECIALIZED AXIS3 PARAMETERS
C	CLOW		REAL VERTICAL DISTANCE BETWEEN CONTOUR AND SURFACE PLOT
C	NLEV		INT  NUMBER OF CONTOUR LEVELS
C			     > 0 NUMBER OF LEVELS IN CVAL
C			     = 0 NO CONTOUR PLOT PRODUCED
C			     < 0 IABS(NLEV) LEVELS COMPUTED
C				 FROM (MAX(A)-MIN(A))/ABS(NLEV-1)
C	CVAL		REAL ARRAY OF CONTOUR LEVELS, IF NLEV<0
C			     CONTOURS USED ARE RETURNED
C	NCH		INT  CONTOUR LABEL OPTION
C			     > 0 LABEL WITH ALPHABETIC CHAR
C			     = 0 NO CONTOUR LABELS
C			     < 0 LABEL WITH CONTOUR VALUE, NCH=# DIGITS
C				 TO RIGHT OF DECIMAL POINT
C	NMIN		INT  MINIMUM NUMBER OF CELLS BEFORE CONTOUR LABELED
C			      < 0 LINE TYPES ARE USED
C			      > 0 LINE TYPES NOT USED
C	ILINE		INT  LINE TYPE/WIDTH ARRAY
C				ILINE(1) = UNDERSIDE OF SURFACE
C				ILINE(2) = INTERCONNECT PATH (RETURN IF USED)
C				ILINE(3) = CONTOUR LINE 1
C				ILINE(4) = CONTOUR LINE 2
C				 ...        (LAST RETURNED IF PATH NOT USED)
C	IPLS		INT  CONTROL FOR LINE PATH INTERCONNECTING SURFACES
C			     AND CONTOUR PLOT
C				IPLS(1) = PLOT CODE POINT 1
C					   = 0 END OF LIST
C					   = 3 MOVE TO, PEN UP
C					   = 2 MOVE TO, PEN DOWN
C				IPLS(2) = X ARRAY INDEX
C				IPLS(3) = Z ARRAY INDEX
C				IPLS(4) = PLOT CODE POINT 2
C				 ...
C FOLLOWING ONLY REQUIRED IF IAXIS IS NOT ZERO
C 	XT,YT,ZT	BYTE STRINGS FOR AXIS TITLES
C 	NXT,NYT,NZT	INT  LENGTH OF AXIS TITLES
C			     IF ZERO THEN THAT AXIS NOT PLOTTED
C 	XASTART,ZASTART	REAL AXIS START VALUES
C 	XAEND,ZAEND	REAL AXIS END VALUES
C FOLLOWING ONLY REQUIRED IF HUNDRED'S DIGIT OF IAXIS=1
C	TX,TY,TZ	REAL TICK MARK PATTERN (SEE AXIS3)
C	TSX,TSY,TSZ	REAL SIZE OF TITLE AND NUMBERS OF AXIS
C			  IF <0 DO NOT AUTO-SCALE AXIS BY (x10^POWER)
C 	FDX,FDY,FDZ	REAL AXIS NUMBER LABEL FORMAT (SEE AXIS3)
C 	AMININ,AMAXIN 	REAL YAXIS SCALING FACTORS (ONLY NEEDED IF IAXIS < 0)
C	ICOL		INTEGER COLOR CONTROL (REQUIRED IF MOD(IFLAG,10)=2)
C				ICOL(1) AXIS LINE COLOR 
C				ICOL(2) AXIS NUMERIC LABEL COLOR 
C				ICOL(3) AXIS LABEL COLOR 
C				ICOL(4) AXIS EXPONENT COLOR 
C				ICOL(5) BACK PLANE COLOR 
C				ICOL(6) PLOT COLOR UPPER SURFACE
C				ICOL(7) PLOT COLOR LOWER SURFACE
C				ICOL(8) CONTOUR LINE 1 COLOR
C				ICOL(9) CONTOUR LINE 2 COLOR
C				 ...     (LAST VALUE RETURNED)
C
	DIMENSION A(INX,INZ),CVAL(1),ILINE(1),IPLS(1),ICOL(1)
	DIMENSION IC(4),AS(2)
	INTEGER XT(1),YT(1),ZT(1)
	LOGICAL REPEAT
	DATA TPI/3.141592654/,REPEAT/.FALSE./
C
C	WORKING SPACE FOR PLT3D
C
	PARAMETER (NW2D0=4000)
	DIMENSION W2(NW2D0)
C
	COMMON /PLT3B/A1,A2,A3,B1,B2,B3,B4
C
	XFORM(I,J,AV)=A1*J+A2*I+A3
	YFORM(I,J,AV)=B1*J+B2*I+B3*AV+B4
C
C	INTIALIZE ROUTINE AND PLOT PACKAGE IF REQUIRED
C
	IF (IFLAG.EQ.0) THEN
		IF (REPEAT) CALL PLOTND		!CLOSE PLOTTER
		REPEAT=.FALSE.
		RETURN
	ENDIF
C
	IF (IAXIS.LT.0) THEN
		AMAX=AMAXIN
		AMIN=AMININ
	ELSE
		AMAX=A(1,1)
		AMIN=A(1,1)
		DO 16IZ=1,NZ		!DETERMINE MAX,MIN ARRAY VALUES
			DO 16 IX=1,NX
				AMAX=AMAX1(AMAX,A(IX,IZ))
				AMIN=AMIN1(AMIN,A(IX,IZ))
16		CONTINUE
	ENDIF
	XLEN=ABS(XH)
	ZLEN=ABS(ZH)
	YLEN=ABS(YH)
	IF (MOD(IABS(IAXIS),10).EQ.2) THEN	!SMOOTH SCALE FACTORS
		AS(1)=AMAX
		AS(2)=AMIN
		CALL SCALE(AS,YLEN,2,1,1,AMIN,DAA)
		AMAX=YLEN*DAA+AMIN
	ENDIF
	YSCALE=1.0
	IF (AMAX-AMIN.NE.0.0) YSCALE=YLEN/(AMAX-AMIN)
C
C	INTIALIZE PLOT PACKAGE
C
	JF=IABS(IFLAG)
	IF (.NOT.REPEAT.AND.JF.LT.10000) THEN
		ILU=JF/100			!GRAPHICS DEVICE CODE
		CALL FRAME(3,ILU,1.5,1.0,1.)	!INTIALIZE LONGLIB
	ELSE
		CALL CTERM(-1)			!PUT TERMINAL IN PLOT MODE
	ENDIF
C
	IAFF=IABS(IAXIS)/100			!AXIS TYPE FLAG
	IAF=MOD(IABS(IAXIS),100)/10		!BACK PANEL, AXES FLAG
	JFF=MOD(JF,100)/10			!SURFACE FLAG
	JF=MOD(JF,10)				!COLOR FLAG
C
C	PLOT MESH SURFACE WITH SIMPLE HIDDEN LINE REMOVAL ROUTINE
C
	NW1D=1
	NW2D=NW2D0-10
	IF (JFF.EQ.0) THEN		! PLT3D
		NW1D=4*MIN(NX,NZ)
		NW2D=NW2D0-NW1D-1
		IF (NW2D.LT.NW1D) GOTO 42
	ENDIF
	NW=NW2D+1
	XP=XLEN/2.*COSD(ALPHA)+2.0
	YP=-(XLEN/2.+YLEN/2.)*SIND(ALPHA)*SIND(BETA+180.)
     $		-YSCALE*AMIN+CLOW
C
C	PLOT LOWER SIDE OF SURFACE
C
	IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(7)),0.,0)	!COLOR SURFACE LOWER
	IF (NMIN.LT.0) CALL NEWPEN(ILINE(1))
	IF (JFF.EQ.0) THEN
	    CALL PLT3D(A,INX,INZ,NX,NZ,W2(NW),NW1D,W2,NW2D,ALPHA,
     $		180.+BETA,-XLEN,XP,ZLEN,YP,YSCALE,0.,IERR)
	ELSE
	    CALL HLT3D(A,INX,INZ,NX,NZ,W2,NW1D,W2,NW2D,ALPHA,
     $		180.+BETA,-XLEN,XP,ZLEN,YP,YSCALE,0.,IERR)
	ENDIF
	IF (IERR.NE.0) GOTO 42
C
C	IF WE ARE SHOWING BACK PANEL, THEN PLOT THE PORTION UNDER SURFACE
C
	IF (IAF.NE.2.AND.IAXIS.NE.0) THEN
		DY=1./YSCALE
		IF (NYT.GT.0) THEN 	!PLOT Y AXIS
			IF (IAFF.EQ.1) THEN
			  N1=INT(TY)	! MAJOR AXIS TICKS
			ELSE
			  N1=IFIX(YLEN+.5)+1
			ENDIF
			IF (IAF.EQ.0.OR.IAF.EQ.4) THEN
			IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(5)),0.,0)
				IF (NMIN.LT.0) CALL NEWPEN(0)
				 VAL2=(N1-1)*DY+AMIN
				 W2(  NW)=XFORM(1,NZ,VAL2)
				 W2(1+NW)=YFORM(1,NZ,VAL2)
				 W2(2+NW)=XFORM(1,1,VAL2)
				 W2(3+NW)=YFORM(1,1,VAL2)
				 W2(4+NW)=XFORM(NX,1,VAL2)
				 W2(5+NW)=YFORM(NX,1,VAL2)
			CALL NXTVU(-2,W2(NW),3,W2,NW2D,IERR)
				 IF (IERR.NE.0) GOTO 42
				DO 338 K=1,N1-1
				 VAL1=(N1-K-1)*DY+AMIN
				 W2(  NW)=XFORM(1,NZ,VAL2)
				 W2(1+NW)=YFORM(1,NZ,VAL2)
				 W2(2+NW)=XFORM(1,NZ,VAL1)
				 W2(3+NW)=YFORM(1,NZ,VAL1)
				 W2(4+NW)=XFORM(1,1,VAL1)
				 W2(5+NW)=YFORM(1,1,VAL1)
				 W2(6+NW)=XFORM(1,1,VAL2)
				 W2(7+NW)=YFORM(1,1,VAL2)
			CALL NXTVU(-2,W2(NW),4,W2,NW2D,IERR)
				 IF (IERR.NE.0) GOTO 42
				 W2(  NW)=XFORM(1,1,VAL1)
				 W2(1+NW)=YFORM(1,1,VAL1)
				 W2(2+NW)=XFORM(NX,1,VAL1)
				 W2(3+NW)=YFORM(NX,1,VAL1)
				 W2(4+NW)=XFORM(NX,1,VAL2)
				 W2(5+NW)=YFORM(NX,1,VAL2)
			CALL NXTVU(-2,W2(NW),3,W2,NW2D,IERR)
				 IF (IERR.NE.0) GOTO 42
338				CONTINUE
			ENDIF
		ENDIF
	ENDIF
C
C	PLOT UPPER SIDE OF SURFACE
C
	IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(6)),0.,0)	!COLOR SURFACE UPPER
	IF (NMIN.LT.0) CALL NEWPEN(0)
	IF (JFF.EQ.0) THEN
	    CALL PLT3D(A,INX,INZ,NX,NZ,W2(NW),NW1D,W2,NW2D,ALPHA,
     $		180.+BETA,XLEN,XP,ZLEN,YP,YSCALE,0.,IERR)
	ELSE
	    CALL HLT3D(A,INX,INZ,NX,NZ,W2,NW1D,W2,NW2D,ALPHA,
     $		180.+BETA,XLEN,XP,ZLEN,YP,YSCALE,0.,IERR)
	ENDIF
	IF (IERR.NE.0) THEN
42		CALL CTERM(1)
		WRITE(*,1) IERR
1	FORMAT(' *** ERROR CALLING PLT3D/HLT3D IN CVAX3DX *** '/
     $		' *** INTERNAL WORKING SPACE EXCEEDED ***',I2)
		RETURN
	ENDIF
	IF (IAF.NE.2.AND.IAXIS.NE.0) THEN
		NADD=0
		IF (JF.EQ.2) THEN
			NADD=100000
			IC(1)=ICOL(1)
			IC(2)=ICOL(2)
			IC(3)=ICOL(3)
			IC(4)=ICOL(4)
		ENDIF
		XP=XFORM(NX,NZ,AMIN)
		YP=YFORM(NX,NZ,AMIN)
		XP1=XFORM(1,NZ,AMIN)
		YP1=YFORM(1,NZ,AMIN)
		ANG=ATAN2(YP-YP1,XP-XP1)*180./TPI
		ZLEN1=SQRT((XP-XP1)**2+(YP-YP1)**2)
		IF (NXT.GT.0) THEN	! PLOT X AXIS
			IF (IAFF.EQ.1) THEN
			  CALL AXIS3(XP1,YP1,XT,-NADD-NXT,ZLEN1,
     $				ANG,XASTART,XAEND,TX,
     $					TSX,FDX,IC)
			ELSE
			  VAL1=INT(XLEN+.5)+1.
			  CALL AXIS3(XP1,YP1,XT,-NADD-NXT,ZLEN1,
     $				ANG,XASTART,XAEND,VAL1,
     $					0.15,5.02,IC)
			ENDIF
		ENDIF
		XP1=XFORM(NX,1,AMIN)
		YP1=YFORM(NX,1,AMIN)
		ANG=ATAN2(YP1-YP,XP1-XP)*180./TPI
		ZLEN1=SQRT((XP1-XP)**2+(YP1-YP)**2)
		IF (NZT.GT.0) THEN		! PLOT Z AXIS
			IF (IAFF.EQ.1) THEN
			  CALL AXIS3(XP,YP,ZT,-NADD-NZT,ZLEN1,
     $				ANG,ZAEND,ZASTART,TZ,
     $					TSZ,FDZ,IC)
			ELSE
			  VAL1=INT(ZLEN+.5)+1.
			  CALL AXIS3(XP,YP,ZT,-NZT,ZLEN1,
     $				ANG,ZAEND,ZASTART,VAL1,
     $					0.15,5.02,IC)
			ENDIF
		ENDIF
C
C	PLOT Y AXISES FOR MESH SURFACE
C
		DY=1./YSCALE
		XP1=XFORM(NX,NZ,AMAX)
		YP1=YFORM(NX,NZ,AMAX)
		ZLEN1=SQRT((XP1-XP)**2+(YP1-YP)**2)
		IF (NYT.GT.0) THEN 	!PLOT Y AXIS
			XP=XFORM(1,NZ,AMIN)	! LEFT SIDE
			YP=YFORM(1,NZ,AMIN)
			XP1=XFORM(NX,1,AMIN)	! RIGHT SIDE
			YP1=YFORM(NX,1,AMIN)
			IF (IAFF.EQ.1) THEN
			  CALL AXIS3(XP,YP,YT,NYT+NADD,ZLEN1,
     $				90.,AMIN,AMAX,TY,
     $					TSY,FDY,IC)
			  CALL AXIS3(XP1,YP1,YT,-NADD-NYT,ZLEN1,
     $					90.,AMIN,AMAX,TY,
     $					TSY,FDY,IC)
			  N1=INT(TY)	! MAJOR AXIS TICKS
			ELSE
			  VAL1=INT(YLEN+.5)+1.
			  CALL AXIS3(XP,YP,YT,NADD+NYT,ZLEN1,
     $				90.,AMIN,AMAX,VAL1,0.15,5.02,IC)
			  CALL AXIS3(XP1,YP1,YT,-NADD-NYT,ZLEN1,
     $				90.,AMIN,AMAX,VAL1,0.15,5.02,IC)
			  N1=IFIX(YLEN+.5)+1
			ENDIF
C
C	PLOT BACK PANEL FOR MESH SURFACE
C
			IF (IAF.EQ.0.OR.IAF.EQ.4) THEN
			IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(5)),0.,0)
				DO 4338 K=1,N1
				 VAL1=(K-1)*DY+AMIN
				 W2(  NW)=XFORM(1,NZ,AMIN)
				 W2(1+NW)=YFORM(1,NZ,AMIN)
				 W2(2+NW)=XFORM(1,NZ,VAL1)
				 W2(3+NW)=YFORM(1,NZ,VAL1)
				 W2(4+NW)=XFORM(1,1,VAL1)
				 W2(5+NW)=YFORM(1,1,VAL1)
				 W2(6+NW)=XFORM(1,1,AMIN)
				 W2(7+NW)=YFORM(1,1,AMIN)
			CALL NXTVU(2,W2(NW),4,W2,NW2D,IERR)
				 IF (IERR.NE.0) GOTO 42
				 W2(  NW)=XFORM(1,1,VAL1)
				 W2(1+NW)=YFORM(1,1,VAL1)
				 W2(2+NW)=XFORM(NX,1,VAL1)
				 W2(3+NW)=YFORM(NX,1,VAL1)
				 W2(4+NW)=XFORM(NX,1,AMIN)
				 W2(5+NW)=YFORM(NX,1,AMIN)
			CALL NXTVU(2,W2(NW),3,W2,NW2D,IERR)
				 IF (IERR.NE.0) GOTO 42
4338				CONTINUE
			ENDIF
		ENDIF
	ENDIF
C
C	DRAW AXES FOR CONTOUR PLOT
C
	IF (NLEV.EQ.0) GOTO 900
	B4=B4-CLOW
	IF (IAF.LT.3.AND.IAXIS.NE.0) THEN
		NADD=0
		IF (JF.EQ.2) THEN
			NADD=100000
			IC(1)=ICOL(1)
			IC(2)=ICOL(2)
			IC(3)=ICOL(3)
			IC(4)=ICOL(4)
		ENDIF
		XP=XFORM(NX,NZ,AMIN)
		YP=YFORM(NX,NZ,AMIN)
		XP1=XFORM(1,NZ,AMIN)
		YP1=YFORM(1,NZ,AMIN)
		ANG=ATAN2(YP-YP1,XP-XP1)*180./TPI
		ZLEN1=SQRT((XP-XP1)**2+(YP-YP1)**2)
		IF (NXT.GT.0) THEN	! PLOT X AXIS
			IF (IAFF.EQ.1) THEN
			  CALL AXIS3(XP1,YP1,XT,-NADD-NXT,ZLEN1,
     $					ANG,XASTART,XAEND,TX,
     $					TSX,FDX,IC)
			ELSE
			  VAL1=INT(XLEN+.5)+1.
			  CALL AXIS3(XP1,YP1,XT,-NADD-NXT,ZLEN1,
     $					ANG,XASTART,XAEND,VAL1,
     $					0.15,5.02,IC)
			ENDIF
		ENDIF
		CALL PLOT(XP1,YP1,3)
		XP1=XFORM(1,1,AMIN)
		YP1=YFORM(1,1,AMIN)
		CALL PLOT(XP1,YP1,2)
		XP1=XFORM(NX,1,AMIN)
		YP1=YFORM(NX,1,AMIN)
		CALL PLOT(XP1,YP1,2)
		ANG=ATAN2(YP1-YP,XP1-XP)*180./TPI
		ZLEN1=SQRT((XP1-XP)**2+(YP1-YP)**2)
		IF (NZT.GT.0) THEN		! PLOT Z AXIS
			IF (IAFF.EQ.1) THEN
			  CALL AXIS3(XP,YP,ZT,-NADD-NZT,ZLEN1,
     $				ANG,ZAEND,ZASTART,TZ,
     $					TSZ,FDZ,IC)
			ELSE
			  VAL1=INT(ZLEN+.5)+1.
			  CALL AXIS3(XP,YP,ZT,-NADD-NZT,ZLEN1,
     $				ANG,ZAEND,ZASTART,VAL1,
     $					0.15,5.02,IC)
			ENDIF
		ENDIF
	ELSE
C
C	BOX CONTOUR PLOT AREA
C
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(1)),0.,0)
		XP=XFORM(NX,NZ,AMIN)
		YP=YFORM(NX,NZ,AMIN)
		CALL PLOT(XP,YP,3)
		XP1=XFORM(1,NZ,AMIN)
		YP1=YFORM(1,NZ,AMIN)
		CALL PLOT(XP1,YP1,2)
		XP1=XFORM(1,1,AMIN)
		YP1=YFORM(1,1,AMIN)
		CALL PLOT(XP1,YP1,2)
		XP1=XFORM(NX,1,AMIN)
		YP1=YFORM(NX,1,AMIN)
		CALL PLOT(XP1,YP1,2)
		CALL PLOT(XP,YP,2)
	ENDIF
C
C	CALL CONTOURING ROUTINE ONCE FOR EACH CONTOUR LEVEL 
C
	B4=B4+B3*AMIN
	AE=1.E20
	ICL=0
	IF (NMIN.LT.0) ICL=2
	NL=IABS(NLEV)
	DO 100 IX=1,NL
		IF (NLEV.LT.1) THEN		! COMPUTE LEVELS
			CVAL(IX)=(IX-1)*(AMAX-AMIN)/(NL-1)+AMIN
		ELSE IF (NLEV.EQ.-1) THEN
			CVAL(IX)=AMIN
		ENDIF
		CV=CVAL(IX)
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(IX+8)),0.,0)	! COLOR
		IF (ICL.NE.0) ILL=ILINE(IX+2)			! LINE TYPE
		CALL GCONTR(A,INX,INZ,NX,NZ,1.,1.,CV,-IX-64,AE,W2,
     $			NCH,-0.12,IABS(NMIN),0,ICL,ILL)
100	CONTINUE
C
C	NOW DRAW INTERCONNECTING PATH LINES
C
	IF (IPLS(1).NE.0) THEN
		IF (JF.EQ.2) CALL PLOT(FLOAT(ICOL(8)),0.,0)	! COLOR
		IF (ICL.NE.0) CALL NEWPEN(ILINE(2))		! LINE TYPE
		ILL=1
		ILAST=1
		JLAST=1
150		CONTINUE
		I=(ILL-1)*3+1
		IF (IPLS(I).LE.0) GOTO 900
		IF (IPLS(I).EQ.2) THEN
			I1=ILAST
			J1=JLAST
			X1=XFORM(I1,J1,A(I1,J1)-AMIN)
			Y1=YFORM(I1,J1,A(I1,J1)-AMIN)+CLOW
			CALL PLOT(X1,Y1,3)
			I1=IPLS(I+1)
			J1=IPLS(I+2)
			X1=XFORM(I1,J1,A(I1,J1)-AMIN)
			Y1=YFORM(I1,J1,A(I1,J1)-AMIN)+CLOW
			CALL PLOT(X1,Y1,2)
			I1=ILAST
			J1=JLAST
			X1=XFORM(I1,J1,0.)
			Y1=YFORM(I1,J1,0.)
			CALL PLOT(X1,Y1,3)
			I1=IPLS(I+1)
			J1=IPLS(I+2)
			X1=XFORM(I1,J1,0.)
			Y1=YFORM(I1,J1,0.)
			CALL PLOT(X1,Y1,2)
		ENDIF
		I1=IPLS(I+1)
		J1=IPLS(I+2)
		X1=XFORM(I1,J1,A(I1,J1)-AMIN)
		Y1=YFORM(I1,J1,A(I1,J1)-AMIN)+CLOW
		CALL PLOT(X1,Y1,3)
		X1=XFORM(I1,J1,0.)
		Y1=YFORM(I1,J1,0.)
		CALL PLOT(X1,Y1,2)
		ILAST=IPLS(I+1)
		JLAST=IPLS(I+2)
		ILL=ILL+1
		GOTO 150
	ENDIF
C
C	FINISH UP
C
900	CALL PLOT(0.,0.,3)			!PEN UP
	IF (IFLAG.GT.0) THEN
		CALL CTERM(2)			!CLEAR TERMINAL
		CALL PLOTND			!CLOSE PLOT OUTPUT
		REPEAT=.FALSE.
	ENDIF
	IF (IFLAG.LT.0) THEN
		REPEAT=.TRUE.
		CALL CTERM(1)			!PUT TERMINAL IN TEXT MODE
	ENDIF
	RETURN
	END
C
C
	SUBROUTINE CNDRAW(X,Y,IL,IH,CVAL,NCH,CS,NMIN,ICL,ICOL,ILINE)
C
C	ROUTINE TO PLOT CONTOUR LINES AND LABELS
C
C	X,Y	(R)	LOCATION OF POINT (IN POINTS)
C	IL	(I)	OPERATION ID FLAG
C			= 1 CONTINUE CONTOUR
C			= 2 START CONTOUR AT BOUNDRY
C			= 3 START CONTOUR NOT AT BOUNDRY
C			= 4 FINISH CONTOUR AT BOUNDRY
C			= 5 FINISH CLOSED CONTOUR (NOT AT BOUNDRY)
C	IH	(I)	CONTOUR NUMBER
C			 IF IH > 0 THEN CONTOUR NUMBER IS LOOKED UP IN
C			 CVAL ARRAY.  IF IH < 0 THEN THE FIRST ELEMENT
C			 OF CVAL IS USED AND CONTOUR IS LABELED WITH IABS(IH)
C	CVAL	(R)	CONTOUR VALUE
C	NCH	(I)	CONTOUR LABELING OPTION
C			< 0 LABEL WITH CONTOUR VALUE (NCH IS NUMBER
C				OF DIGITS TO THE LEFT OF DECIMAL
C			= 0 NO LABELING OF CONTOURS
C			> 0 LABEL WITH LETTERS
C	CS	(R)	LABEL SIZE
C			< 0 PLT3D LOCATION SCALING USED WITH Z VALUE
C			    OF ZERO
C			> 0 NORMAL SCALING USED
C	NMIN	(I)	MINIMUM NUMBER OF SEGMENTS TO BE LABELED
C	ICL	(I)	COLOR AND LINE TYPE FLAG
C			= 0 NO COLOR OR LINE TYPE USED
C			= 1 COLOR USED ONLY
C			= 2 LINE TYPE USED ONLY
C			= 3 COLOR AND LINE TYPE USED
C	ICOL	(I)	LINE COLOR FOR EACH CONTOUR
C	ILINE	(I)	LINE TYPE FOR EACH CONTOUR
C
	DIMENSION ILINE(1),ICOL(1)
	INTEGER WORK(1)
C
	COMMON /PLT3B/A1,A2,A3,B1,B2,B3,B4
C
	IPEN=2
	IF (IL.EQ.2.OR.IL.EQ.3) THEN
		IHN=MAX(1,IH)
		IF (ICL.GT.1) CALL NEWPEN(ILINE(IHN))	! LINE TYPE
		IF (MOD(ICL,2).EQ.1) CALL PLOT(FLOAT(ICOL(IHN)),0.,0)! PEN COLOR
		IPEN=3
		NCNT=0
	ENDIF
	NCNT=NCNT+1
	X1=X
	Y1=Y
	IF (CS.LT.0.0) THEN
		X1=A1*(Y+1.0)+A2*(X+1.0)+A3
		Y1=B1*(Y+1.0)+B2*(X+1.0)+B4
	ENDIF
	CALL PLOT(X1,Y1,IPEN)
	IF (IL.LT.4) GOTO 30
	IF (NCH.EQ.0) GOTO 30
	IF (NCNT.LT.NMIN) GOTO 30
	IF (ICL.GT.1) CALL NEWPEN(0)			! RESET LINE TYPE
	IF (NCH.LT.0) GOTO 5
	IF (IABS(IH).LT.26) GOTO 10
5	CALL NUMBER(X1,Y1-0.05,ABS(CS),CVAL,0.,IABS(NCH),-1)
	GOTO 30
10	CALL SYMBOL(X1,Y1-0.05,ABS(CS),64+IABS(IH),0.,1,-1)
30	RETURN
	END
C
C
	SUBROUTINE GCONTR(Z,NDX,NDY,NX,NY,XL,YL,CV,NC,ZMAX,BITMAP,NCH,CS,
     $		NMIN,ICL,ICOL,ILINE)
C
C	DRAWS CONTOURS THROUGH EQUAL VALUES OF AN ORTHOGONAL ARRAY
C
C	Z	(R)	ARRAY FOR WHICH CONTOURS DRAWN
C			DIMENSIONED Z(NDX,NDY) (X,Y)
C	NDX,NDY	(I)	DIMENSION OF Z
C	NX,NY	(I)	NUMBER OF X,Y POINTS IN Z TO USE (NX<=NRZ)
C	XL,YL	(R)	SCALE FACTORS FOR X,Y POINTS (INCHES/POINT)
C	CV	(R)	ARRAY OF CONTOUR VALUES TO BE DRAWN DIM DV(NCV)
C			CV < 0 CAUSES PLT3D PLOTTING PARAMETERS TO BE USED
C			FOR OUTPUTTING CONTOUR VALUES
C	NC	(I)	NUMBER OF CONTOURS IN CV (SET TO NCV)
C			 IF NCV < 0 THEN ONLY ONE CONTOUR IS ASSUMED.  IT
C			 WILL BE LABELED AS CONTOUR IABS(NC) AND NCV WILL BE
C			 EQUAL TO 1 WHEREVER IT APPEARS.
C	ZMAX	(R)	MAXIMUM VALUE OF Z FOR CONSIDERATION.  A POINT
C			Z(I,J) WHIC EXCEEDS ZMAX WILL CAUSE THIS POINT AND
C			THE LINE SEGMENTS WHICH RADIATE OUT FROM IT TO NOT
C			BE CONSIDERED WHEN CONTOURING.
C	BITMAP	(I*4)	WORK SPACE DIMENSIONED AT LEAST
C			(2*NX*NY*NCV+NBPW-1)/NBPW WHERE NBPW IS THE
C			NUMBER OF USEFUL BITS IN A WORD (NBPW=31)
C	NCH	(I)	CONTOUR LABELING OPTION
C			< 0 LABEL WITH CONTOUR VALUE (NCH IS NUMBER
C				OF DIGITS TO THE LEFT OF DECIMAL)
C			= 0 NO LABELING OF CONTOURS
C			> 0 LABEL WITH LETTERS
C	CS	(R)	SIZE OF LABELS
C			> 0 NORMAL LINEAR OPERATION
C			< 0 PLT3D TRANSFORMATION USED (XL,YL SHOULD BE 1.0)
C	NMIN	(I)	MINIMUM NUMBER OF SEGMENTS CROSSED FOR CONTOUR
C			TO BE LABELED
C	ICL	(I)	COLOR AND LINE TYPE FLAG
C			= 0 NO COLOR OR LINE TYPE USED
C			= 1 COLOR USED ONLY
C			= 2 LINE TYPE USED ONLY
C			= 3 COLOR AND LINE TYPE USED
C	ICOL	(I)	LINE COLOR FOR EACH CONTOUR
C	ILINE	(I)	LINE TYPE FOR EACH CONTOUR
C
	REAL Z(NDX,NDY),CV(1)
	INTEGER ICOL(1),ILINE(1)
	INTEGER*4 BITMAP(1),IBIT(32)
	INTEGER L1(4),L2(4),IJ(2)	! LIMITS USED IN SPIRAL SEARCH
	INTEGER I1(2),I2(2),I3(6)	! NEIGHBOR COMPUTATIONS
	REAL XINT(4)			! INTERSECTIONS OF CONTOURS WITH CELL EDGES
	REAL XY(2)			! COORDINATES FOR PLOTTING
	EQUIVALENCE (L2(1),IMAX),(L2(2),JMAX),(L2(3),IMIN),(L2(4),JMIN)
	EQUIVALENCE (IJ(1),I),(IJ(2),J),(XY(1),X),(XY(2),Y)
C
	DATA L1(3)/-1/,L1(4)/-1/,I1/1,0/,I2/1,-1/,I3/1,0,0,1,1,0/
	DATA NBPW/31/
C
C	IGET RETURNS THE VALUE OF THE NTH BIT OF BITMAP
C
	IGET(N)=MOD(BITMAP((N-1)/NBPW+1)/2**(NBPW-MOD(N-1,NBPW)-1),2)
C
C	COMPUTE BITS ARRAY
C
	NCV=1
	IF (NC.GT.0) NCV=NC
	DO 2 I=1,NBPW
		IBIT(I)=2**(I-1)
2	CONTINUE
C
	L1(1)=NX
	L1(2)=NY
	DMAX=ZMAX
C
C	SET INITIAL POSITION
C
	X=NX
	Y=NY
	ICUR=X
	JCUR=Y
C
C	CLEAR BITMAP
C	FILLS BITMAP WITH N ZEROS.  NBPW IS THE NUMBER OF USEFUL BITS IN
C	A WORD (USUALLY 1 LESS THAN THE NUMBER OF BITS)
C
	N=2*NX*NY*NCV
	LOOP=N/NBPW
	NBLW=MOD(N,NBPW)
	IF (LOOP.EQ.0) GOTO 5
	DO 1 I=1,LOOP
		BITMAP(I)=0
1	CONTINUE
5	IF(NBLW.NE.0) BITMAP(LOOP+1)=0
C
C	SEARCH ALONG RECTANGULAR SPIRAL PATH FOR A LINE SEGMENT THAT
C	1. END POINTS ARE NOT EXCLUDED
C	2. NO MARK HAS BEEN RECORDED FOR THE SEGMENT
C	3. VALUES OF Z AT END OF SEGMENT ARE SUTCH THAT ONE Z IS LESS THAN
C	   THE CURRENT CONTOUR VALUE, AND THE OTHER IS GREATER THAN OR EQUAL
C	   TO THE CURRENT CONTOUR VALUE
C
	IBKEY=0
10	I=ICUR
	J=JCUR
20	IMAX=I
	IMIN=-I
	JMAX=J
	JMIN=-J
	IDIR=0
C
C	DIRECTION CODE:
C	 0: +I   1:+J   2:-I   3:-J
C
30	NXIDIR=IDIR+1
	K=NXIDIR
	IF (NXIDIR.GT.3) NXIDIR=0
40	I=IABS(I)
	J=IABS(J)
	IF (Z(I,J).GT.DMAX) GOTO 140
	L=1
C L=1 MEANS HORIZONTAL LINE, L=2 MEANS VERTICAL LINE
50	IF (IJ(L).GE.L1(L)) GOTO 130
	II=I+I1(L)
	JJ=J+I1(3-L)
	IF (Z(II,JJ).GT.DMAX) GOTO 130
C
	ASSIGN 100 TO JUMP
C
C THE NEXT 15 STATEMENTS DETECT BOUNDRIES
C
60	IX=1
	IF (IJ(3-L).EQ.1) GOTO 80
	II=I-I1(3-L)
	JJ=J-I1(L)
	IF (Z(II,JJ).GT.DMAX) GOTO 70
	II=I+I2(L)
	JJ=J+I2(3-L)
	IF (Z(II,JJ).LT.DMAX) IX=0
70	IF (IJ(3-L).GE.L1(3-L)) GOTO 90
80	II=I+I1(3-L)
	JJ=J+I1(L)
	IF (Z(II,JJ).GT.DMAX) GOTO 90
	IF (Z(I+1,J+1).LT.DMAX) GOTO JUMP 	! (100 OR 280)
90	IX=IX+2
	GOTO JUMP 				! (100 OR 280)
100	IF (IX.EQ.3) GOTO 130
	IF (IX+IBKEY.EQ.0) GOTO 130
C DETERMINE WHETHER LINE SEGMENT IS CROSSED BY CONTOUR
	II=I+I1(L)
	JJ=J+I1(3-L)
	Z1=Z(I,J)
	Z2=Z(II,JJ)
	DO 120 ICV=1,NCV
CC		IF (IGET(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L).NE.0) GOTO 120
		IF (IGET(2*(NX*(NY*(ICV-1)+J-1)+I-1)+L).NE.0) GOTO 120
		IF (CV(ICV).LE.MIN(Z1,Z2)) GOTO 110
		IF (CV(ICV).LE.MAX(Z1,Z2)) GOTO 190
110		CONTINUE
C
C	MARK BIT MAP
C
CC		CALL MARK1(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L,NBPW)
		LL=2*(NX*(NY*(ICV-1)+J-1)+I-1)+L
		NWORD=(LL-1)/NBPW+1
		NBIT=MOD(LL-1,NBPW)
		IIJ=IBIT(NBPW-NBIT)
		BITMAP(NWORD)=BITMAP(NWORD)+IIJ*(1-MOD(BITMAP(NWORD)/IIJ,2))
120	CONTINUE
130	L=L+1
	IF (L.LE.2) GOTO 50
140	L=MOD(IDIR,2)+1
	IJ(L)=ISIGN(IJ(L),L1(K))
C
C	LINES FROM Z(I,J) TO Z(I+1,J) AND Z(I,J+1) ARE NOT SATISFACTORY,
C	SO CONTINUE SPIRAL
C
150	IF (IJ(L).GE.L1(K)) GOTO 170
	IJ(L)=IJ(L)+1
	IF (IJ(L).GT.L2(K)) GOTO 160
	GOTO 40
160	L2(K)=IJ(L)
	IDIR=NXIDIR
	GOTO 30
170	IF (IDIR.EQ.NXIDIR) GOTO 180
	NXIDIR=NXIDIR+1
	IJ(L)=L1(K)
	K=NXIDIR
	L=3-L
	IJ(L)=L2(K)
	IF (NXIDIR.GT.3) NXIDIR=0
	GOTO 150
180	IF (IBKEY.NE.0) RETURN
	IBKEY=1
	GOTO 10
C
C AN ACCEPTABLE LINE SEGMENT HAS BEEN FOUND.  FOLLOW CONTOUR UNTIL IT 
C HITS THE BOUNDRY OR CLOSES
C
190	IEDGE=L
	CVAL=CV(ICV)
	IF (IX.NE.1) IEDGE=IEDGE+2
	IFLAG=2+IBKEY
	XINT(IEDGE)=(CVAL-Z1)/(Z2-Z1)
200	XY(L)=FLOAT(IJ(L))+XINT(IEDGE)
	XY(3-L)=FLOAT(IJ(3-L))
C
C	MARK BIT MAP
C
CC	CALL MARK1(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L,NBPW)
	LL=2*(NX*(NY*(ICV-1)+J-1)+I-1)+L
	NWORD=(LL-1)/NBPW+1
	NBIT=MOD(LL-1,NBPW)
	IIJ=IBIT(NBPW-NBIT)
	BITMAP(NWORD)=BITMAP(NWORD)+IIJ*(1-MOD(BITMAP(NWORD)/IIJ,2))
C
	XX=(X-1.)*XL
	YY=(Y-1.)*YL
	CALL CNDRAW(XX,YY,IFLAG,NC,CVAL,NCH,CS,NMIN,ICL,ICOL,ILINE)
	IF (IFLAG.LT.4) GOTO 210
	ICUR=I
	JCUR=J
	GOTO 20
C
C CONTINUE CONTOUR.  EDGES ARE NUMBERED CLOCKWISE WITH BOTTOM EDGE BEING
C EDGE NUMBER ONE.
C
210	NI=1
	IF (IEDGE.LT.3) GOTO 220
	I=I-I3(IEDGE)
	J=J-I3(IEDGE+2)
220	DO 250 K=1,4
		IF (K.EQ.IEDGE) GOTO 250
		II=I+I3(K)
		JJ=J+I3(K+1)
		Z1=Z(II,JJ)
		II=I+I3(K+1)
		JJ=J+I3(K+2)
		Z2=Z(II,JJ)
		IF (CVAL.LE.AMIN1(Z1,Z2)) GOTO 250
		IF (CVAL.GT.AMAX1(Z1,Z2)) GOTO 250
		IF (K.EQ.1) GOTO 230
		IF (K.NE.4) GOTO 240
230		ZZ=Z1
		Z1=Z2
		Z2=ZZ
240		XINT(K)=(CVAL-Z1)/(Z2-Z1)
		NI=NI+1
		KS=K
250	CONTINUE
	IF (NI.EQ.2) GOTO 260
C
C CONTOUR CROSSES ALL FOUR EDGES OF THE CELL.  CHOOSE LINES TOP-TO-LEFT
C AND BOTTOM TO RIGHT IF THE INTERPOLATION POINT ON THE TOP EDGE IS LESS THAN
C THE INTERPOLATION POINT ON THE BOTTOM EDGE.  OTHERWISE CHOSE THE OTHER
C PAIR.  METHOD PRODUCES SAME RESULTS FOR REVERSED AXES.  CONTOUR MAY CLOSE
C AT ANY EDGE, BUT MUST NOT CROSS ITSELF IN ANY CELL.
C
	KS=5-IEDGE
	IF (XINT(3).LT.XINT(1)) GOTO 260
	KS=3-IEDGE
	IF (KS.LE.0) KS=KS+4
C
260	L=KS
	IFLAG=1
	ASSIGN 280 TO JUMP
	IF (KS.LT.3) GOTO 270
	I=I+I3(KS)
	J=J+I3(KS+2)
	L=KS-2
270	IF (IGET(2*(NX*(NY*(ICV-1)+J-1)+I-1)+L).EQ.0) GOTO 60
	IFLAG=5
	GOTO 290
280	IF (IX.NE.0) IFLAG=4
290	IEDGE=KS+2
	IF (IEDGE.GT.4) IEDGE=IEDGE-4
	XINT(IEDGE)=XINT(KS)
	GOTO 200
	END
